---
title: "\\Large Risk of a feedback loop between climatic warming and human mobility"
subtitle: "Supplementary Information"
author: Nick Obradovich and Iyad Rahwan
date: ""
output: 
  pdf_document:
    fig_caption: true
    toc: false
    keep_tex: true
fontsize: 11pt
header-includes:
    - \usepackage{caption}
    - \usepackage{bbm}
    - \makeatletter\renewcommand*{\fps@figure}{h}\makeatother
---

```{r opts, echo=FALSE, message=FALSE, warning=FALSE, include=FALSE}
knitr::opts_chunk$set(fig.path='./', warning=F, message=F, fig.retina=NULL, cache=T, cache.lazy=F, autodep=T, echo=F, eval=T)
```

```{r packages}
# packages----
library(bit64) # long integers
library(ggthemes) #themes for plotting
library(sandwich) #robust se estimator
library(lmtest) #lm functions
library(ordinal) #ologit funcitons
library(memisc) #import stata data
library(ggmap) #make nice plots with ggplot mapping
library(ggplot2) #nice plots
library(Hmisc) #multipurpose package
library(RColorBrewer) #make nice colors
library(grid) #plotting package
library(MASS) #multipurpose package
library(gplots) #nice plotting
library(lfe) #run linear fixed effects models
library(stargazer) #output nice tables
library(lubridate) #deal with dates
library(zoo) #deal with dates
library(reshape2) #powerful reshape package
library(sp) #general spatial class package
library(gstat) #geostat package
library(spacetime) #store data in proper spatial/temporal class
library(raster) #convert to raster and vice versa
library(maptools) #plot and modify shapefiles
library(rgeos) #for use with maptools for shapefiles
library(rgdal) #for reading in shapefiles
library(RSAGA) #needed for spatial downscaling
library(RCurl) #also needed for spatial downscaling
library(dplyr) #for data frame manipulation functions
library(plyr) #for additional manipulation functions
library(data.table) #for data.table functions
library(foreign) #for reading in .dta files
library(parallel) #for parallel processing functions
library(foreach) #for plyr parallel
library(doMC) #run multicore processes
library(caTools) #functions for fast running mean and sd
library(ncdf4) #tools for netCDF packages
library(stringr) #string tools
library(knitr) #knitting
library(pander) #pandering
library(rasterVis) # plotting rasters
library(GGally) # for correlation plots
library(geosphere) # spherical distance
library(Rcpp) # for C++
library(RcppArmadillo) # for C++ (and Conley errors)
library(matrixStats) # matrix calculations
library(ggExtra) # for marginal histograms
library(gridExtra) # for multiplots
library(viridis) # for colors
library(Cairo) # pdf device
library(gtable) # for gtable filter
library(devtools)
#devtools::install_github("wilkelab/cowplot")
library(cowplot)
library(captioner)
library(splines)
library(survey) # svycontrast

# set cores for parallel processing----
registerDoMC(20)
```

```{r functions}
# functions for grabbing p-values and coefficients from summaries of felm models----
coefs <- function(reg, name) {
    if (round(reg$coefficients[rownames(reg$coefficients)==name, 1], 3) == 0) {
        format(round(reg$coefficients[rownames(reg$coefficients)==name, 1], 4), scientific=FALSE) 
    }
    else round(reg$coefficients[rownames(reg$coefficients)==name, 1], 3)
}
pval <- function(reg, name) {
  if (round(reg$coefficients[rownames(reg$coefficients)==name, 4], 3) >=0.001) 
  {round(reg$coefficients[rownames(reg$coefficients)==name, 4], 3)} 
  else paste("< 0.001")
}

# captioner setting----
fig_num <- captioner(prefix="Fig.")
tab_num <- captioner(prefix="Table")
eq_num <- captioner(prefix="Equation")
sfig_num <- captioner(prefix="Fig. S", auto_space=FALSE)
stab_num <- captioner(prefix="Table S", auto_space=FALSE)
seq_num <- captioner(prefix="Equation S", auto_space=FALSE)

# binned plot----
## convert felm result into a data frame with rhs | coef | se | xmin | xmax | xmid | ci.l | ci.h,
## replacing Inf and -Inf with in xmin and xmax on the edge bins with the graphically spaced bound
## felm.est: result from felm() call (lm will probably work too) with RHS factor variable like cuttmax or bins
## plotvar: string of variable to plot: e.g. 'cuttemp'
## breaks: bin size
## omit: omitted bin
bin.plot <- function(felm.est, plotvar, breaks, omit, panel="panel", ylabel="y-label", xlabel="x=label",
                     y.axis.title = element_text(size=25, angle=90, margin = margin(t = 0, r = 10, b = 0, l = 0), face="bold"),
                     ylimit = c(min(dt$ci.l) - (max(dt$ci.h) - min(dt$ci.l))/10,max(dt$ci.h)),
                     y.axis.text = element_text(size=13, angle=0, face="bold"),
                     y.axis.line = element_line(),
                     y.axis.ticks = element_line(),
                     errorline="gray60", errorfill="transparent", linecolor="#FB8C00", pointfill="#FB8C00", pointcolor="white") {

    dt <- as.data.frame(summary(felm.est)$coefficients[, 1:2]) # extract from felm summary (use df to keep coefficient names)
    dt <- dt[grepl(plotvar, rownames(dt)), ] # extract variable name (e.g., tmax.cut, ppt.cut, ...), limit to matches
    dt$rhs <- str_split_fixed(rownames(dt), "   ", n=2)[,1]
    dt <- as.data.table(dt)
    dt[, rhs := gsub(paste0(plotvar), "", rhs)]
    dt[, rhs := gsub("\\(|]", "", rhs)]
    dt[, rhs := gsub("neg", "-", rhs)]
    dt[, rhs := gsub("pos", "", rhs)]
    dt[grepl(",", rhs), rhs := tstrsplit(rhs, ',')[2]] # if a factor variable, grab second column (upper limit)
    dt[, rhs := as.numeric(rhs)]
    dt <- dt[!is.na(rhs)]
    dt[, xmin := rhs - breaks]
    dt[, xmax := rhs]
    dt[, rhs := NULL]
    setnames(dt, c("coef", "se", "xmin", "xmax"))
    dt[, xmin.lead := data.table::shift(xmin, type="lead")]
    dt[, xmax.lag := data.table::shift(xmax, type="lag")]
    dt[xmin==-Inf, xmin := xmin.lead - breaks]
    dt[xmax==-Inf, xmax := xmin.lead]
    dt[xmin==Inf, xmin := xmax.lag]
    dt[xmax==Inf, xmax := xmin + breaks]
    dt[, c('xmin.lead','xmax.lag') := NULL]
    # identify omitted category, add 0.
    dt <- rbind(dt, data.table(coef=0, se=0, xmin=omit[1], xmax=omit[2]))
    dt[, ':='(xmid=(xmin+xmax)/2, ci.l=coef - se*1.96, ci.h=coef + se*1.96)] # add some stuff for plotting
    dt[, range := paste(xmin, xmax, sep="-")] # create range variable for x labels
    dt[xmax == min(xmax), range := paste0("<=",min(xmax))] # set bottom of range
    dt[xmin == max(xmin), range := paste0(">",min(xmin))] # set top of range

    plot <- ggplot() +
        geom_hline(aes(yintercept=0), linetype=2, color='gray60') +
        geom_ribbon(data=dt, aes(x=xmid, ymin=ci.l, ymax=ci.h), fill=errorfill, alpha=0.1) +
        geom_line(data=dt, aes(x=xmid, y=ci.l), linetype="dotted", color=errorline, size=1.5, alpha=1) +
        geom_line(data=dt, aes(x=xmid, y=ci.h), linetype="dotted", color=errorline, size=1.5, alpha=1) +
        geom_line(data=dt, aes(x=xmid, y=coef), color=linecolor, size=1.75, alpha=1) +
        geom_point(data=dt, aes(x=xmid, y=coef), size=3, fill=pointfill, color=pointcolor, pch=21, alpha=1) +
        annotate("text", label=panel, x = min(dt$xmid) + (max(dt$xmid) - min(dt$xmid))/10, y = min(dt$ci.l) - (max(dt$ci.h) - min(dt$ci.l))/10, size=9) + # plot panel annotation
        ylab(ylabel) +
        xlab(xlabel) +
        coord_cartesian(xlim=c(min(dt$xmid), max(dt$xmid)), ylim=ylimit, expand=TRUE) +
        scale_x_continuous(breaks = dt$xmid, labels = dt$range) +
        theme(plot.title = element_text(face="bold", size=40),
              axis.title.x = element_text(size=20, angle=0, face="bold"),
              axis.title.y = y.axis.title,
              axis.text.y = y.axis.text,
              axis.text.x = element_text(size=13, face="bold"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.x = element_line(),
              axis.line.y = y.axis.line,
              axis.ticks.y = y.axis.ticks,
              legend.title=element_blank()
        )
    return(plot)
}

# spline functions----
## vmt
vmt_estimate_spline <- function(data, dv, iv, lhs) {
  ## keep spline knots with the fitted model results for plotting
  spl_knots <- wtd.quantile(iv[!is.na(dv)], weights=dv[!is.na(dv)], probs = c(0.25, 0.5, 0.75), na.rm=T) # get quantiles of iv weighted by dv (cannot have missing values)
  spl_boundary_knots <- range(iv, na.rm=T) # get range of independent variable
  xs <- as.data.table(ns(iv, knots = spl_knots, Boundary.knots = spl_boundary_knots)) # get cubic splines of iv
  setnames(xs, paste0("ns_", 1:ncol(xs))) # set names of splines
  data <- cbind(data, xs)
  ## set up felm model
  rhs <- paste(c(names(xs), "trange", "prcp.days", "cloud", "humid", "wind"), collapse = " + ")
  fe <- paste(c("state:month", "state:year"), collapse = " + ")
  fmla <- as.formula(sprintf("%s ~ %s | %s | 0 | state", lhs, rhs, fe))
  fit_felm <- felm(fmla, data=data)
  pred_tmax <- seq(min(iv, na.rm=T), max(iv, na.rm=T), length.out = 200) # prediction iv
  xp <- as.data.table(ns(pred_tmax, knots = spl_knots, Boundary.knots = spl_boundary_knots)) 
  result <- foreach::foreach (i = 1:nrow(xp), .combine = rbind) %do% {
    this.x <- pred_tmax[i]
    this.row <- c(unlist(xp[i, ]), 0, 0, 0, 0, 0) # include zeros for meteorological controls
    temp.dt <- as.data.table(svycontrast(fit_felm, this.row))
    temp.dt[, x := this.x]
    setnames(temp.dt, c("coef", "se", "x"))
    return(temp.dt)
  }
  return(result)
} # modified function from patrick baylis

## upt
upt_estimate_spline <- function(data, dv, iv, lhs) {
  ## keep spline knots with the fitted model results for plotting
  spl_knots <- wtd.quantile(iv[!is.na(dv)], weights=dv[!is.na(dv)], probs = c(0.25, 0.5, 0.75), na.rm=T) # get quantiles of iv weighted by dv (cannot have missing values)
  spl_boundary_knots <- range(iv, na.rm=T) # get range of independent variable
  xs <- as.data.table(ns(iv, knots = spl_knots, Boundary.knots = spl_boundary_knots)) # get cubic splines of iv
  setnames(xs, paste0("ns_", 1:ncol(xs))) # set names of splines
  data <- cbind(data, xs)
  ## set up felm model
  rhs <- paste(c(names(xs), "trange", "prcp.days", "cloud", "humid", "wind"), collapse = " + ")
  fe <- paste(c("uza:month", "uza:year"), collapse = " + ")
  fmla <- as.formula(sprintf("%s ~ %s | %s | 0 | uza", lhs, rhs, fe))
  fit_felm <- felm(fmla, data=data)
  pred_tmax <- seq(min(iv, na.rm=T), max(iv, na.rm=T), length.out = 200) # prediction iv
  xp <- as.data.table(ns(pred_tmax, knots = spl_knots, Boundary.knots = spl_boundary_knots)) 
  result <- foreach::foreach (i = 1:nrow(xp), .combine = rbind) %do% {
    this.x <- pred_tmax[i]
    this.row <- c(unlist(xp[i, ]), 0, 0, 0, 0, 0) # include zeros for meteorological controls
    temp.dt <- as.data.table(svycontrast(fit_felm, this.row))
    temp.dt[, x := this.x]
    setnames(temp.dt, c("coef", "se", "x"))
    return(temp.dt)
  }
  return(result)
} # modified function from patrick baylis
```

```{r data}
load('./replication.RData')
```

```{r models}
# vehicle miles traveled main bootstrap----
vmt.main <- readRDS("./vmt/main_bootstrap.rds")

# public transit trips main bootstrap----
upt.main <- readRDS("./upt/main_bootstrap.rds")
```

\captionsetup[table]{labelformat=empty}
\captionsetup[figure]{labelformat=empty}

<!-- \listoffigures -->
<!-- \listoftables -->
\newpage

<!-- figures-->

```{r fig.hetero, eval=FALSE}
# urban vmt and rural vmt bootstrap----
## urban
# boot.results <- foreach::foreach (b = 1:1000, .combine = rbind) %dopar% {
#                     bstrap_sample <- vmt[.(sample(unique(state), replace = T))] # block bootstrap sample on state
#                     a <- vmt_estimate_spline(data=bstrap_sample, dv=bstrap_sample$log_urban_vmt, iv=bstrap_sample$tmax, lhs="log_urban_vmt")[, b := b] # estimate spline
#                     return(a)
#                   }
# main.result <- vmt_estimate_spline(data=vmt, dv=vmt$log_urban_vmt, iv=vmt$tmax, lhs="log_urban_vmt")[, b:= 0] # main
# boot.results <- rbind(boot.results, main.result) # combine main and bootstrapped results
# boot.results[, bmax := max(coef, na.rm=T), by = b] # get max of that bootstrap's prediction
# boot.results[, coef_norm := coef - bmax] # normalize so max=0
# saveRDS(boot.results, "./vmt/urban_bootstrap.rds")
# rm(main.result, boot.results)
vmt.urban <- readRDS("./vmt/urban_bootstrap.rds")

## rural
# boot.results <- foreach::foreach (b = 1:1000, .combine = rbind) %dopar% {
#                     bstrap_sample <- vmt[.(sample(unique(state), replace = T))] # block bootstrap sample on state
#                     a <- vmt_estimate_spline(data=bstrap_sample, dv=bstrap_sample$log_rural_vmt, iv=bstrap_sample$tmax, lhs="log_rural_vmt")[, b := b] # estimate spline
#                     return(a)
#                   }
# main.result <- vmt_estimate_spline(data=vmt, dv=vmt$log_rural_vmt, iv=vmt$tmax, lhs="log_rural_vmt")[, b:= 0] # main
# boot.results <- rbind(boot.results, main.result) # combine main and bootstrapped results
# boot.results[, bmax := max(coef, na.rm=T), by = b] # get max of that bootstrap's prediction
# boot.results[, coef_norm := coef - bmax] # normalize so max=0
# saveRDS(boot.results, "./vmt/rural_bootstrap.rds")
# rm(main.result, boot.results)
vmt.rural <- readRDS("./vmt/rural_bootstrap.rds")

# split vmt by median average tmax bootstrap----
## hot
# boot.results <- foreach::foreach (b = 1:1000, .combine = rbind) %dopar% {
#                     bstrap_sample <- vmt[.(sample(unique(state), replace = T))][warm==TRUE] # block bootstrap sample on state
#                     a <- vmt_estimate_spline(data=bstrap_sample, dv=bstrap_sample$log_vmt, iv=bstrap_sample$tmax, lhs="log_vmt")[, b := b] # estimate spline
#                     return(a)
#                   }
# main.result <- vmt_estimate_spline(data=vmt[warm==TRUE], dv=vmt[warm==TRUE, log_vmt], iv=vmt[warm==TRUE, tmax], lhs="log_vmt")[, b:= 0] # main
# boot.results <- rbind(boot.results, main.result) # combine main and bootstrapped results
# boot.results[, bmax := max(coef, na.rm=T), by = b] # get max of that bootstrap's prediction
# boot.results[, coef_norm := coef - bmax] # normalize so max=0
# saveRDS(boot.results, "./vmt/hot_bootstrap.rds")
# rm(main.result, boot.results)
vmt.hot <- readRDS("./vmt/hot_bootstrap.rds")

## cold
# boot.results <- foreach::foreach (b = 1:1000, .combine = rbind) %dopar% {
#                     bstrap_sample <- vmt[.(sample(unique(state), replace = T))][warm==FALSE] # block bootstrap sample on state
#                     a <- vmt_estimate_spline(data=bstrap_sample, dv=bstrap_sample$log_vmt, iv=bstrap_sample$tmax, lhs="log_vmt")[, b := b] # estimate spline
#                     return(a)
#                   }
# main.result <- vmt_estimate_spline(data=vmt[warm==FALSE], dv=vmt[warm==FALSE, log_vmt], iv=vmt[warm==FALSE, tmax], lhs="log_vmt")[, b:= 0] # main
# boot.results <- rbind(boot.results, main.result) # combine main and bootstrapped results
# boot.results[, bmax := max(coef, na.rm=T), by = b] # get max of that bootstrap's prediction
# boot.results[, coef_norm := coef - bmax] # normalize so max=0
# saveRDS(boot.results, "./vmt/cold_bootstrap.rds")
# rm(main.result, boot.results)
vmt.cold <- readRDS("./vmt/cold_bootstrap.rds")

# split vmt by median of vmt per capita bootstrap----
## high
# boot.results <- foreach::foreach (b = 1:1000, .combine = rbind) %dopar% {
#                     bstrap_sample <- vmt[.(sample(unique(state), replace = T))][high==TRUE] # block bootstrap sample on state
#                     a <- vmt_estimate_spline(data=bstrap_sample, dv=bstrap_sample$log_vmt, iv=bstrap_sample$tmax, lhs="log_vmt")[, b := b] # estimate spline
#                     return(a)
#                   }
# main.result <- vmt_estimate_spline(data=vmt[high==TRUE], dv=vmt[high==TRUE, log_vmt], iv=vmt[high==TRUE, tmax], lhs="log_vmt")[, b:= 0] # main
# boot.results <- rbind(boot.results, main.result) # combine main and bootstrapped results
# boot.results[, bmax := max(coef, na.rm=T), by = b] # get max of that bootstrap's prediction
# boot.results[, coef_norm := coef - bmax] # normalize so max=0
# saveRDS(boot.results, "./vmt/highvmt_bootstrap.rds")
# rm(main.result, boot.results)
vmt.high <- readRDS("./vmt/highvmt_bootstrap.rds")

## low
# boot.results <- foreach::foreach (b = 1:1000, .combine = rbind) %dopar% {
#                     bstrap_sample <- vmt[.(sample(unique(state), replace = T))][high==FALSE] # block bootstrap sample on state
#                     a <- vmt_estimate_spline(data=bstrap_sample, dv=bstrap_sample$log_vmt, iv=bstrap_sample$tmax, lhs="log_vmt")[, b := b] # estimate spline
#                     return(a)
#                   }
# main.result <- vmt_estimate_spline(data=vmt[high==FALSE], dv=vmt[high==FALSE, log_vmt], iv=vmt[high==FALSE, tmax], lhs="log_vmt")[, b:= 0] # main
# boot.results <- rbind(boot.results, main.result) # combine main and bootstrapped results
# boot.results[, bmax := max(coef, na.rm=T), by = b] # get max of that bootstrap's prediction
# boot.results[, coef_norm := coef - bmax] # normalize so max=0
# saveRDS(boot.results, "./vmt/lowvmt_bootstrap.rds")
# rm(main.result, boot.results)
vmt.low <- readRDS("./vmt/lowvmt_bootstrap.rds")

# split upt by rail vs nonrail bootstrap----
## rail
# boot.results <- foreach::foreach (b = 1:1000, .combine = rbind) %dopar% {
#                     k <- 0
#                     while(k!=1000){
#                           bstrap_sample = try(upt[.(sample(unique(uza), replace = T))])
#                           if (class(bstrap_sample)=="try-error") {
#                              k <- k+1
#                              print(k)
#                              } else break
#                           } # block bootstrap sample on city, don't oversample (code handles errors)
#                     a <- upt_estimate_spline(data=bstrap_sample, dv=bstrap_sample$log_rail_trips, iv=bstrap_sample$tmax, lhs="log_rail_trips")[, b := b] # estimate spline
#                     return(a)
#                   }
# main.result <- upt_estimate_spline(data=upt, dv=upt$log_rail_trips, iv=upt$tmax, lhs="log_rail_trips")[, b:= 0] # main
# boot.results <- rbind(boot.results, main.result) # combine main and bootstrapped results
# boot.results[, bmax := max(coef, na.rm=T), by = b] # get max of that bootstrap's prediction
# boot.results[, coef_norm := coef - bmax] # normalize so max=0
# saveRDS(boot.results, "./upt/rail_bootstrap.rds")
# rm(main.result, boot.results)
upt.rail <- readRDS("./upt/rail_bootstrap.rds")

## nonrail
# boot.results <- foreach::foreach (b = 1:1000, .combine = rbind) %dopar% {
#                     k <- 0
#                     while(k!=1000){
#                           bstrap_sample = try(upt[.(sample(unique(uza), replace = T))])
#                           if (class(bstrap_sample)=="try-error") {
#                              k <- k+1
#                              print(k)
#                              } else break
#                           } # block bootstrap sample on city, don't oversample (code handles errors)
#                     a <- upt_estimate_spline(data=bstrap_sample, dv=bstrap_sample$log_nonrail_trips, iv=bstrap_sample$tmax, lhs="log_nonrail_trips")[, b := b] # estimate spline
#                     return(a)
#                   }
# main.result <- upt_estimate_spline(data=upt, dv=upt$log_nonrail_trips, iv=upt$tmax, lhs="log_nonrail_trips")[, b:= 0] # main
# boot.results <- rbind(boot.results, main.result) # combine main and bootstrapped results
# boot.results[, bmax := max(coef, na.rm=T), by = b] # get max of that bootstrap's prediction
# boot.results[, coef_norm := coef - bmax] # normalize so max=0
# saveRDS(boot.results, "./upt/nonrail_bootstrap.rds")
# rm(main.result, boot.results)
upt.nonrail <- readRDS("./upt/nonrail_bootstrap.rds")

# split upt by median average tmax bootstrap----
## hot
# boot.results <- foreach::foreach (b = 1:1000, .combine = rbind) %dopar% {
#                     k <- 0
#                     while(k!=1000){
#                           bstrap_sample = try(upt[.(sample(unique(uza), replace = T))][warm==TRUE])
#                           if (class(bstrap_sample)=="try-error") {
#                              k <- k+1
#                              print(k)
#                              } else break
#                           } # block bootstrap sample on city, don't oversample (code handles errors)
#                     a <- upt_estimate_spline(data=bstrap_sample, dv=bstrap_sample$log_trips, iv=bstrap_sample$tmax, lhs="log_trips")[, b := b] # estimate spline
#                     return(a)
#                   }
# main.result <- upt_estimate_spline(data=upt[warm==TRUE], dv=upt[warm==TRUE, log_trips], iv=upt[warm==TRUE, tmax], lhs="log_trips")[, b:= 0] # main
# boot.results <- rbind(boot.results, main.result) # combine main and bootstrapped results
# boot.results[, bmax := max(coef, na.rm=T), by = b] # get max of that bootstrap's prediction
# boot.results[, coef_norm := coef - bmax] # normalize so max=0
# saveRDS(boot.results, "./upt/hot_bootstrap.rds")
# rm(main.result, boot.results)
upt.hot <- readRDS("./upt/hot_bootstrap.rds")

## cold
# boot.results <- foreach::foreach (b = 1:1000, .combine = rbind) %dopar% {
#                     k <- 0
#                     while(k!=1000){
#                           bstrap_sample = try(upt[.(sample(unique(uza), replace = T))][warm==FALSE])
#                           if (class(bstrap_sample)=="try-error") {
#                              k <- k+1
#                              print(k)
#                              } else break
#                           } # block bootstrap sample on city, don't oversample (code handles errors)
#                     a <- upt_estimate_spline(data=bstrap_sample, dv=bstrap_sample$log_trips, iv=bstrap_sample$tmax, lhs="log_trips")[, b := b] # estimate spline
#                     return(a)
#                   }
# main.result <- upt_estimate_spline(data=upt[warm==FALSE], dv=upt[warm==FALSE, log_trips], iv=upt[warm==FALSE, tmax], lhs="log_trips")[, b:= 0] # main
# boot.results <- rbind(boot.results, main.result) # combine main and bootstrapped results
# boot.results[, bmax := max(coef, na.rm=T), by = b] # get max of that bootstrap's prediction
# boot.results[, coef_norm := coef - bmax] # normalize so max=0
# saveRDS(boot.results, "./upt/cold_bootstrap.rds")
# rm(main.result, boot.results)
upt.cold <- readRDS("./upt/cold_bootstrap.rds")

# split upt by median of transit trips per capita 2016 bootstrap----
## high
# boot.results <- foreach::foreach (b = 1:1000, .combine = rbind) %dopar% {
#                     k <- 0
#                     while(k!=1000){
#                           bstrap_sample = try(upt[.(sample(unique(uza), replace = T))][high==TRUE])
#                           if (class(bstrap_sample)=="try-error") {
#                              k <- k+1
#                              print(k)
#                              } else break
#                           } # block bootstrap sample on city, don't oversample (code handles errors)
#                     a <- upt_estimate_spline(data=bstrap_sample, dv=bstrap_sample$log_trips, iv=bstrap_sample$tmax, lhs="log_trips")[, b := b] # estimate spline
#                     return(a)
#                   }
# main.result <- upt_estimate_spline(data=upt[high==TRUE], dv=upt[high==TRUE, log_trips], iv=upt[high==TRUE, tmax], lhs="log_trips")[, b:= 0] # main
# boot.results <- rbind(boot.results, main.result) # combine main and bootstrapped results
# boot.results[, bmax := max(coef, na.rm=T), by = b] # get max of that bootstrap's prediction
# boot.results[, coef_norm := coef - bmax] # normalize so max=0
# saveRDS(boot.results, "./upt/hightrips_bootstrap.rds")
# rm(main.result, boot.results)
upt.high <- readRDS("./upt/hightrips_bootstrap.rds")

## low
# boot.results <- foreach::foreach (b = 1:1000, .combine = rbind) %dopar% {
#                     k <- 0
#                     while(k!=1000){
#                           bstrap_sample = try(upt[.(sample(unique(uza), replace = T))][high==FALSE])
#                           if (class(bstrap_sample)=="try-error") {
#                              k <- k+1
#                              print(k)
#                              } else break
#                           } # block bootstrap sample on city, don't oversample (code handles errors)
#                     a <- upt_estimate_spline(data=bstrap_sample, dv=bstrap_sample$log_trips, iv=bstrap_sample$tmax, lhs="log_trips")[, b := b] # estimate spline
#                     return(a)
#                   }
# main.result <- upt_estimate_spline(data=upt[high==FALSE], dv=upt[high==FALSE, log_trips], iv=upt[high==FALSE, tmax], lhs="log_trips")[, b:= 0] # main
# boot.results <- rbind(boot.results, main.result) # combine main and bootstrapped results
# boot.results[, bmax := max(coef, na.rm=T), by = b] # get max of that bootstrap's prediction
# boot.results[, coef_norm := coef - bmax] # normalize so max=0
# saveRDS(boot.results, "./upt/lowtrips_bootstrap.rds")
# rm(main.result, boot.results)
upt.low <- readRDS("./upt/lowtrips_bootstrap.rds")

# vehicle miles traveled urban/rural----
dum <- data.frame(point = sort(rep(letters[1:2], 10)), x = runif(20), y = runif(20))
vmt.urbanrural.plot <- ggplot() + 
                        geom_hline(aes(yintercept=0), linetype=2, color='gray80') +
                        geom_line(data=vmt.urban[b!=0], aes(x=x, y=coef_norm, group=b), color="gray30", alpha=0.02) +
                        geom_line(data=vmt.rural[b!=0], aes(x=x, y=coef_norm, group=b), color="#FB8C00", alpha=0.02) + 
                        geom_line(data=vmt.urban[b==0,], aes(x=x, y=coef_norm), color="gray30", size=3, alpha=0.9) +
                        geom_line(data=vmt.rural[b==0,], aes(x=x, y=coef_norm), color="#FB8C00", size=3, alpha=0.9) +
                        coord_cartesian(ylim=c(-0.25, 0.01)) +
                        scale_x_continuous(breaks = seq(-10,40,10)) +
                        ylab("Log of Vehicle Miles Traveled") +
                        xlab("") +
                        geom_line(data = dum, aes(color = point, x = x, y = y), alpha = 0, size=3) + 
                        scale_color_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                           limits=c("gray30", "#FB8C00"), 
                           values=c("gray30", "#FB8C00"),
                           labels=c("Urban", "Rural")) +
                        theme(plot.title = element_text(face="bold", size=40),
                              axis.title.x = element_text(size=25, face="bold"),
                              axis.title.y = element_blank(),
                              axis.text.y = element_blank(),
                              axis.text.x = element_text(size=13, face="bold"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.background = element_blank(),
                              axis.line.x = element_line(),
                              axis.line.y = element_blank(),
                              axis.ticks.y = element_blank(),
                              legend.title=element_blank(),
                              legend.position = c(.4, 0.2),
                              legend.text=element_text(size=16, face="bold"),
                              legend.key = element_rect(colour = "black", fill="white")
                        )
hist <- axis_canvas(vmt.urbanrural.plot, axis = "x") + 
        geom_histogram(data = vmt, aes(x = tmax), bins=30, fill="#FB8C00", color='gray60', alpha=0.5, size=0.1)
vmt.urbanrural.plot <- insert_xaxis_grob(vmt.urbanrural.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
vmt.urbanrural.plot <- ggdraw(vmt.urbanrural.plot)

# vehicle miles traveled by climate----
dum <- data.frame(point = sort(rep(letters[1:2], 10)), x = runif(20), y = runif(20))
vmt.climate.plot <- ggplot() + 
                        geom_hline(aes(yintercept=0), linetype=2, color='gray80') +
                        geom_line(data=vmt.cold[b!=0], aes(x=x, y=coef_norm, group=b), color="gray50", alpha=0.02) +
                        geom_line(data=vmt.hot[b!=0], aes(x=x, y=coef_norm, group=b), color="#c66e00", alpha=0.02) + 
                        geom_line(data=vmt.cold[b==0,], aes(x=x, y=coef_norm), color="gray60", size=3, alpha=0.9) +
                        geom_line(data=vmt.hot[b==0,], aes(x=x, y=coef_norm), color="#c66e00", size=3, alpha=0.9) +
                        coord_cartesian(ylim=c(-0.25, 0.01)) +
                        scale_x_continuous(breaks = seq(-10,40,10)) +
                        ylab("Log of Vehicle Miles Traveled") +
                        xlab("") +
                        geom_line(data = dum, aes(color = point, x = x, y = y), alpha = 0, size=3) + 
                        scale_color_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                           limits=c("gray50", "#c66e00"), 
                           values=c("gray50", "#c66e00"),
                           labels=c("Cooler", "Warmer")) +
                        theme(plot.title = element_text(face="bold", size=40),
                              axis.title.x = element_text(size=25, face="bold"),
                              axis.title.y = element_blank(),
                              axis.text.y = element_blank(),
                              axis.text.x = element_text(size=13, face="bold"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.background = element_blank(),
                              axis.line.x = element_line(),
                              axis.line.y = element_blank(),
                              axis.ticks.y = element_blank(),
                              legend.title=element_blank(),
                              legend.position = c(.4, 0.2),
                              legend.text=element_text(size=16, face="bold"),
                              legend.key = element_rect(colour = "black", fill="white")
                        )
hist <- axis_canvas(vmt.climate.plot, axis = "x") + 
        geom_histogram(data = vmt, aes(x = tmax), bins=30, fill="#c66e00", color='gray60', alpha=0.5, size=0.1)
vmt.climate.plot <- insert_xaxis_grob(vmt.climate.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
vmt.climate.plot <- ggdraw(vmt.climate.plot)

# vehicle miles traveled by vmt per capita----
dum <- data.frame(point = sort(rep(letters[1:2], 10)), x = runif(20), y = runif(20))
vmt.percap.plot <- ggplot() + 
                        geom_hline(aes(yintercept=0), linetype=2, color='gray80') +
                        geom_line(data=vmt.low[b!=0], aes(x=x, y=coef_norm, group=b), color="gray70", alpha=0.02) + 
                        geom_line(data=vmt.high[b!=0], aes(x=x, y=coef_norm, group=b), color="#ffce91", alpha=0.02) + 
                        geom_line(data=vmt.low[b==0,], aes(x=x, y=coef_norm), color="gray70", size=3, alpha=0.9) +
                        geom_line(data=vmt.high[b==0,], aes(x=x, y=coef_norm), color="#ffce91", size=3, alpha=0.9) +
                        coord_cartesian(ylim=c(-0.25, 0.01)) +
                        scale_x_continuous(breaks = seq(-10,40,10)) +
                        ylab("Log of Vehicle Miles Traveled") +
                        xlab("") +
                        geom_line(data = dum, aes(color = point, x = x, y = y), alpha = 0, size=3) + 
                        scale_color_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                           limits=c("gray70", "#ffce91"), 
                           values=c("gray70", "#ffce91"),
                           labels=c("Lower VMT", "Higher VMT")) +
                        theme(plot.title = element_text(face="bold", size=40),
                              axis.title.x = element_text(size=25, face="bold"),
                              axis.title.y = element_blank(),
                              axis.text.y = element_blank(),
                              axis.text.x = element_text(size=13, face="bold"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.background = element_blank(),
                              axis.line.x = element_line(),
                              axis.line.y = element_blank(),
                              axis.ticks.y = element_blank(),
                              legend.title=element_blank(),
                              legend.position = c(.4, 0.2),
                              legend.text=element_text(size=16, face="bold"),
                              legend.key = element_rect(colour = "black", fill="white")
                        )
hist <- axis_canvas(vmt.percap.plot, axis = "x") + 
        geom_histogram(data = vmt, aes(x = tmax), bins=30, fill="#ffce91", color='gray60', alpha=0.5, size=0.1)
vmt.percap.plot <- insert_xaxis_grob(vmt.percap.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
vmt.percap.plot <- ggdraw(vmt.percap.plot)

# justy.vmt----
justy.vmt <- ggplot() + 
              geom_hline(aes(yintercept=0), linetype=2, color='gray80') +
              geom_line(data=vmt.high[b!=0], aes(x=x, y=coef_norm, group=b), color="#FB8C00", alpha=0.02) + 
              geom_line(data=vmt.low[b!=0], aes(x=x, y=coef_norm, group=b), color="gray60", alpha=0.02) + 
              geom_line(data=vmt.high[b==0,], aes(x=x, y=coef_norm), color="#FB8C00", size=3, alpha=0.9) +
              geom_line(data=vmt.low[b==0,], aes(x=x, y=coef_norm), color="gray60", size=3, alpha=0.9) +
              coord_cartesian(ylim=c(-0.25, 0.01)) +
              scale_x_continuous(breaks = seq(-10,40,10)) +
              ylab("Log(VMT)") +
              xlab("") +
              geom_line(data = dum, aes(color = point, x = x, y = y), alpha = 0, size=3) + 
              scale_color_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                 limits=c("#FB8C00", "gray60"), 
                 values=c("#FB8C00", "gray60"),
                 labels=c("Lower VMT", "Higher VMT")) +
              theme(plot.title = element_text(face="bold", size=40),
                    axis.title.x = element_text(size=25, face="bold"),
                    axis.title.y = element_text(size=25, face="bold"),
                    axis.text.y = element_text(size=13, face="bold"),
                    axis.text.x = element_text(size=13, face="bold"),
                    panel.grid.major = element_blank(),
                    panel.grid.minor = element_blank(),
                    panel.background = element_blank(),
                    axis.line.x = element_line(),
                    axis.line.y = element_line(),
                    axis.ticks.y = element_line(),
                    legend.title=element_blank(),
                    legend.position = c(.4, 0.2),
                    legend.text=element_text(size=16, face="bold"),
                    legend.key = element_rect(colour = "black", fill="white")
              )
justy.vmt <- gtable_filter(ggplotGrob(justy.vmt), 'axis-l|ylab', trim=F)
justy.vmt <- ggdraw(justy.vmt)

# pubic transit trips rail/nonrail----
dum <- data.frame(point = sort(rep(letters[1:2], 10)), x = runif(20), y = runif(20))
upt.railnonrail.plot <- ggplot() + 
                        geom_hline(aes(yintercept=0), linetype=2, color='gray80') +
                        geom_line(data=upt.nonrail[b!=0], aes(x=x, y=coef_norm, group=b), color="gray30", alpha=0.02) +
                        geom_line(data=upt.rail[b!=0], aes(x=x, y=coef_norm, group=b), color="#19A1FC", alpha=0.02) + 
                        geom_line(data=upt.nonrail[b==0,], aes(x=x, y=coef_norm), color="gray30", size=3, alpha=0.9) +
                        geom_line(data=upt.rail[b==0,], aes(x=x, y=coef_norm), color="#19A1FC", size=3, alpha=0.9) +
                        coord_cartesian(ylim=c(-0.5, 0.01)) +
                        scale_x_continuous(breaks = seq(-10,40,10)) +
                        ylab("Log of Public Transit Trips") +
                        xlab("") +
                        geom_line(data = dum, aes(color = point, x = x, y = y), alpha = 0, size=3) + 
                        scale_color_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                           limits=c("#19A1FC", "gray30"), 
                           values=c("#19A1FC", "gray30"),
                           labels=c("Rail", "Non-Rail")) +
                        theme(plot.title = element_text(face="bold", size=40),
                              axis.title.x = element_text(size=25, face="bold"),
                              axis.title.y = element_blank(),
                              axis.text.y = element_blank(),
                              axis.text.x = element_text(size=13, face="bold"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.background = element_blank(),
                              axis.line.x = element_line(),
                              axis.line.y = element_blank(),
                              axis.ticks.y = element_blank(),
                              legend.title=element_blank(),
                              legend.position = c(.4, 0.2),
                              legend.text=element_text(size=16, face="bold"),
                              legend.key = element_rect(colour = "black", fill="white")
                        )
hist <- axis_canvas(upt.railnonrail.plot, axis = "x") + 
        geom_histogram(data = upt, aes(x = tmax), bins=30, fill="#19A1FC", color='gray60', alpha=0.5, size=0.1)
upt.railnonrail.plot <- insert_xaxis_grob(upt.railnonrail.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
upt.railnonrail.plot <- ggdraw(upt.railnonrail.plot)

# public transit trips by climate----
dum <- data.frame(point = sort(rep(letters[1:2], 10)), x = runif(20), y = runif(20))
upt.climate.plot <- ggplot() + 
                        geom_hline(aes(yintercept=0), linetype=2, color='gray80') +
                        geom_line(data=upt.hot[b!=0], aes(x=x, y=coef_norm, group=b), color="gray50", alpha=0.02) + 
                        geom_line(data=upt.cold[b!=0], aes(x=x, y=coef_norm, group=b), color="#0366a8", alpha=0.02) + 
                        geom_line(data=upt.hot[b==0,], aes(x=x, y=coef_norm), color="gray50", size=3, alpha=0.9) +
                        geom_line(data=upt.cold[b==0,], aes(x=x, y=coef_norm), color="#0366a8", size=3, alpha=0.9) +
                        coord_cartesian(ylim=c(-0.5, 0.01)) +
                        scale_x_continuous(breaks = seq(-10,40,10)) +
                        ylab("Log of Public Transit Trips") +
                        xlab("") +
                        geom_line(data = dum, aes(color = point, x = x, y = y), alpha = 0, size=3) + 
                        scale_color_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                           limits=c("#0366a8", "gray50"), 
                           values=c("#0366a8", "gray50"),
                           labels=c("Cooler", "Warmer")) +
                        theme(plot.title = element_text(face="bold", size=40),
                              axis.title.x = element_text(size=25, face="bold"),
                              axis.title.y = element_blank(),
                              axis.text.y = element_blank(),
                              axis.text.x = element_text(size=13, face="bold"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.background = element_blank(),
                              axis.line.x = element_line(),
                              axis.line.y = element_blank(),
                              axis.ticks.y = element_blank(),
                              legend.title=element_blank(),
                              legend.position = c(.4, 0.2),
                              legend.text=element_text(size=16, face="bold"),
                              legend.key = element_rect(colour = "black", fill="white")
                        )
hist <- axis_canvas(upt.climate.plot, axis = "x") + 
        geom_histogram(data = upt, aes(x = tmax), bins=30, fill="#0366a8", color='gray60', alpha=0.5, size=0.1)
upt.climate.plot <- insert_xaxis_grob(upt.climate.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
upt.climate.plot <- ggdraw(upt.climate.plot)

# public transit trips by trips per capita 2016----
dum <- data.frame(point = sort(rep(letters[1:2], 10)), x = runif(20), y = runif(20))
upt.percap.plot <- ggplot() + 
                        geom_hline(aes(yintercept=0), linetype=2, color='gray80') +
                        geom_line(data=upt.high[b!=0], aes(x=x, y=coef_norm, group=b), color="gray70", alpha=0.02) + 
                        geom_line(data=upt.low[b!=0], aes(x=x, y=coef_norm, group=b), color="#6ac1fc", alpha=0.02) + 
                        geom_line(data=upt.high[b==0,], aes(x=x, y=coef_norm), color="gray70", size=3, alpha=0.9) +
                        geom_line(data=upt.low[b==0,], aes(x=x, y=coef_norm), color="#6ac1fc", size=3, alpha=0.9) +
                        coord_cartesian(ylim=c(-0.5, 0.01)) +
                        scale_x_continuous(breaks = seq(-10,40,10)) +
                        ylab("Log of Public Transit Trips") +
                        xlab("") +
                        geom_line(data = dum, aes(color = point, x = x, y = y), alpha = 0, size=3) + 
                        scale_color_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                           limits=c("#6ac1fc", "gray70"), 
                           values=c("#6ac1fc", "gray70"),
                           labels=c("Lower Trips", "Higher Trips")) +
                        theme(plot.title = element_text(face="bold", size=40),
                              axis.title.x = element_text(size=25, face="bold"),
                              axis.title.y = element_blank(),
                              axis.text.y = element_blank(),
                              axis.text.x = element_text(size=13, face="bold"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.background = element_blank(),
                              axis.line.x = element_line(),
                              axis.line.y = element_blank(),
                              axis.ticks.y = element_blank(),
                              legend.title=element_blank(),
                              legend.position = c(.4, 0.2),
                              legend.text=element_text(size=16, face="bold"),
                              legend.key = element_rect(colour = "black", fill="white")
                        )
hist <- axis_canvas(upt.percap.plot, axis = "x") + 
        geom_histogram(data = upt, aes(x = tmax), bins=30, fill="#6ac1fc", color='gray60', alpha=0.5, size=0.1)
upt.percap.plot <- insert_xaxis_grob(upt.percap.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
upt.percap.plot <- ggdraw(upt.percap.plot)

# justy.upt----
justy.upt <- ggplot() + 
              geom_hline(aes(yintercept=0), linetype=2, color='gray80') +
              geom_line(data=upt.high[b!=0], aes(x=x, y=coef_norm, group=b), color="gray60", alpha=0.02) + 
              geom_line(data=upt.low[b!=0], aes(x=x, y=coef_norm, group=b), color="#19A1FC", alpha=0.02) + 
              geom_line(data=upt.high[b==0,], aes(x=x, y=coef_norm), color="gray60", size=3, alpha=0.9) +
              geom_line(data=upt.low[b==0,], aes(x=x, y=coef_norm), color="#19A1FC", size=3, alpha=0.9) +
              coord_cartesian(ylim=c(-0.5, 0.01)) +
              scale_x_continuous(breaks = seq(-10,40,10)) +
              ylab("Log(PTT)") +
              xlab("") +
              geom_line(data = dum, aes(color = point, x = x, y = y), alpha = 0, size=3) + 
              scale_color_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                 limits=c("#19A1FC", "gray60"), 
                 values=c("#19A1FC", "gray60"),
                 labels=c("Lower Trips", "Higher Trips")) +
              theme(plot.title = element_text(face="bold", size=40),
                    axis.title.x = element_text(size=25, face="bold"),
                    axis.title.y = element_text(size=25, face="bold"),
                    axis.text.y = element_text(size=13, face="bold"),
                    axis.text.x = element_text(size=13, face="bold"),
                    panel.grid.major = element_blank(),
                    panel.grid.minor = element_blank(),
                    panel.background = element_blank(),
                    axis.line.x = element_line(),
                    axis.line.y = element_line(),
                    axis.ticks.y = element_line(),
                    legend.title=element_blank(),
                    legend.position = c(.4, 0.2),
                    legend.text=element_text(size=16, face="bold"),
                    legend.key = element_rect(colour = "black", fill="white")
              )
justy.upt <- gtable_filter(ggplotGrob(justy.upt), 'axis-l|ylab', trim=F)
justy.upt <- ggdraw(justy.upt)

# plot----
png(filename = './heterogeneous_effects.png', width=14, height=7, units = "in", type="cairo", res=300)
grid.arrange(
  arrangeGrob(justy.vmt, vmt.urbanrural.plot, vmt.climate.plot, vmt.percap.plot, ncol=4, widths = c(0.07, 0.31, 0.31, 0.31)),
  arrangeGrob(justy.upt, upt.railnonrail.plot, upt.climate.plot, upt.percap.plot, ncol=4, widths = c(0.07, 0.31, 0.31, 0.31)),
  nrow=2
)
grid.text("a", x=unit(0.08, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("b", x=unit(0.4, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("c", x=unit(0.71, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("d", x=unit(0.08, "npc"), y=unit(0.495, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("e", x=unit(0.4, "npc"), y=unit(0.495, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("f", x=unit(0.71, "npc"), y=unit(0.495, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("Average Maximum Monthly Temperature in \u2103", x=unit(0.55, "npc"), y=unit(0.035, "npc"), gp=gpar(fontsize=25,font=2))
dev.off()
```

```{r hetero, out.width="100%", fig.cap = paste0(sfig_num("hetero", display="cite"),": Plot of heterogeneous effects. This figure plots the heterogeneous tests for each of VMT and PTT. Panel (a) plots the effects of temperature on VMT in urban centers as compared to these effects on VMT in rural areas. Panel (b) breaks out the sample of states along median average maximum temperatures and conducts the splined regressions for each half of the sample. Panel (c) calculates average VMT per capita in our sample and splits the sample along the median of this metric. Across panels (a)-(c) we observe little evidence of heterogeneous effects; temperatures appear to have very similar effects on VMT across context. Panel (e) breaks out PTT by those that occur via rail versus those that occur via non-rail means (primarily bus travel). Rail travel responds more acutely to extreme temperatures, though non-rail travel also decreases in cold temperatures. Panel (f) splits the sample by median average maximum temperatures and indicates that warmer states have more acute responses to colder temperatures than do colder states. Panel (g) splits the sample by median PTT per capita and indicates that temperature impacts states similarly across baseline transit usage. In each panel, error regions are given by conducting and plotting 1,000 cluster bootstrapped spline regressions for each sample.")}
knitr::include_graphics("./heterogeneous_effects.png")
```

```{r}
hetero <- sfig_num(name = "hetero", caption = "")
```

```{r fig.robustness, eval=FALSE}
# daybin count models----
## vmt
vmt.bin.main <- felm(log_vmt ~ tmaxneginf + tmax5 + tmax10 + tmax15 + tmax20 + tmax25 + tmax30 
                  + trange + prcp.days + cloud + humid + wind 
                  | state:month + state:year 
                  | 0 
                  | state
                  , data=vmt)
vmt.main.bin.plot <- bin.plot(felm.est = vmt.bin.main,
                          plotvar = "tmax",
                          breaks = 5,
                          omit = c(30,35),
                          ylimit = c(-0.006, 0.0005),
                          xlabel = "Max. Temp. in \u2103",
                          ylabel = "Log(VMT)",
                          panel = "")
hist <- axis_canvas(vmt.main.bin.plot, axis = "x") + 
  geom_histogram(data = vmt, aes(x = tmax), bins=30, fill="#FB8C00", color='gray60', alpha=0.5, size=0.1)
vmt.main.bin.plot <- insert_xaxis_grob(vmt.main.bin.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
vmt.main.bin.plot <- ggdraw(vmt.main.bin.plot)

## upt
upt.bin.main <- felm(log_trips ~ tmaxneginf + tmax5 + tmax10 + tmax15 + tmax20 + tmax25 + tmax30 
                  + trange + prcp.days + cloud + humid + wind 
                  | uza:month + uza:year 
                  | 0 
                  | uza
                  , data=upt)
upt.main.bin.plot <- bin.plot(felm.est = upt.bin.main,
                          linecolor = "#19A1FC",
                          pointfill = "#19A1FC",
                          plotvar = "tmax",
                          breaks = 5,
                          omit = c(30,35),
                          ylimit = c(-0.008, 0.0005),
                          xlabel = "Max. Temp. in \u2103",
                          ylabel = "Log(PTT)",
                          panel = "")
hist <- axis_canvas(upt.main.bin.plot, axis = "x") + 
  geom_histogram(data = upt, aes(x = tmax), bins=30, fill="#19A1FC", color='gray60', alpha=0.5, size=0.1)
upt.main.bin.plot <- insert_xaxis_grob(upt.main.bin.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
upt.main.bin.plot <- ggdraw(upt.main.bin.plot)

# monthly tmax binned models----
## vmt
vmt[, cuttmax := cut(tmax, breaks=c(-Inf,seq(0, 30, by=5),Inf))] # create factor variable for regresison
vmt[, cuttmax := relevel(cuttmax, ref="(30, Inf]")] # set reference category

vmt.mon.bin <- felm(log_vmt ~ cuttmax
                  + trange + prcp.days + cloud + humid + wind 
                  | state:month + state:year 
                  | 0 
                  | state
                  , data=vmt)
vmt.mon.bin.plot <- bin.plot(felm.est = vmt.mon.bin,
                          plotvar = "tmax",
                          breaks = 5,
                          omit = c(30,35),
                          ylimit = c(-0.1, 0.005),
                          xlabel = "Max. Temp. in \u2103",
                          ylabel = "Log(VMT)",
                          panel = "")
hist <- axis_canvas(vmt.mon.bin.plot, axis = "x") + 
  geom_histogram(data = vmt, aes(x = tmax), bins=30, fill="#FB8C00", color='gray60', alpha=0.5, size=0.1)
vmt.mon.bin.plot <- insert_xaxis_grob(vmt.mon.bin.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
vmt.mon.bin.plot <- ggdraw(vmt.mon.bin.plot)

## upt
upt[, cuttmax := cut(tmax, breaks=c(-Inf,seq(0, 30, by=5),Inf))] # create factor variable for regresison
upt[, cuttmax := relevel(cuttmax, ref="(30, Inf]")] # set reference category

upt.mon.bin <- felm(log_trips ~ cuttmax
                  + trange + prcp.days + cloud + humid + wind 
                  | uza:month + uza:year 
                  | 0 
                  | uza
                  , data=upt)
upt.mon.bin.plot <- bin.plot(felm.est = upt.mon.bin,
                          linecolor = "#19A1FC",
                          pointfill = "#19A1FC",
                          plotvar = "tmax",
                          breaks = 5,
                          omit = c(30,35),
                          ylimit = c(-0.15, 0.005),
                          xlabel = "Max. Temp. in \u2103",
                          ylabel = "Log(PTT)",
                          panel = "")
hist <- axis_canvas(upt.mon.bin.plot, axis = "x") + 
  geom_histogram(data = upt, aes(x = tmax), bins=30, fill="#19A1FC", color='gray60', alpha=0.5, size=0.1)
upt.mon.bin.plot <- insert_xaxis_grob(upt.mon.bin.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
upt.mon.bin.plot <- ggdraw(upt.mon.bin.plot)

# quadratic monthly tmax models----
## vmt
vmt[, tmax2 := tmax^2]
vmt.quad <- felm(log_vmt ~ tmax + tmax2
                  + trange + prcp.days + cloud + humid + wind 
                  | state:month + state:year 
                  | 0 
                  | state
                  , data=vmt)

t <- seq(from=-10, to=40, by=0.001)
d.t <- vmt.quad$coef[1]*t + vmt.quad$coef[2]*t^2
d.t <- d.t - max(d.t)
vmt.t.data <- as.data.table(cbind(t, d.t))

vmt.t.plot <- ggplot(vmt.t.data) + 
    ylab("Log(VMT)") +
    xlab("Max. Temp. in \u2103") +
    geom_hline(yintercept = 0, linetype = 2, color="gray60") +
    geom_line(aes(x=t, y=d.t), color="#FB8C00", size=1.5) + 
    scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
    coord_cartesian(xlim=c(-10,40), ylim = c(-0.2, 0.01)) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=20, face="bold"),
          axis.title.y = element_text(size=25, face="bold"),
          axis.text.y = element_text(size=13, face="bold"),
          axis.text.x = element_text(size=13, face="bold"),
          axis.ticks.y = element_line(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_line(),
          legend.title=element_blank(),
          legend.position = "none"
    )
hist <- axis_canvas(vmt.t.plot, axis = "x") + 
        geom_histogram(data = vmt, aes(x = tmax), bins=30, fill='#FB8C00', color='gray60', alpha=0.5, size=0.1)
vmt.t.plot <- insert_xaxis_grob(vmt.t.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
vmt.t.plot <- ggdraw(vmt.t.plot)

##  margn effect temp
vmt.quad.vcv <- vmt.quad$clustervcv # get cluster robust variance covariance matrices

# get coefficients
vmt.quad.tcoef <- vmt.quad$coef[1]
vmt.quad.t2coef <- vmt.quad$coef[2]

# calculate marginal effects
marg <- data.table(tmax = seq(from=-10, to=40, by=0.001))
marg[, vmt.quad.delta := vmt.quad.tcoef + 2*vmt.quad.t2coef*tmax]

# compute variance of marginal effect
marg[, vmt.quad.se := sqrt(vmt.quad.vcv["tmax","tmax"] + 4*(tmax^2)*vmt.quad.vcv["tmax2","tmax2"] + 4*tmax*vmt.quad.vcv[1,2])]

# 95% confidence bounds
marg[, vmt.quad.upr := vmt.quad.delta + 1.96*vmt.quad.se]
marg[, vmt.quad.lwr := vmt.quad.delta - 1.96*vmt.quad.se]

vmt.quad.plot <- ggplot(data=marg) + 
                ylab("Marginal Effect of +1\u2103\non Log(VMT)") +
                xlab("Max. Temp. in \u2103") +
                geom_hline(aes(yintercept=0), linetype="dashed", color="gray60") +
                geom_ribbon(aes(x=tmax, ymin=vmt.quad.lwr, ymax=vmt.quad.upr), fill="gray60", alpha=0.25) + 
                geom_line(aes(x=tmax, y=vmt.quad.delta), color="#FB8C00", size=1.5) +
                scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
                coord_cartesian(ylim=c(-0.002, 0.01)) +
                theme(plot.title = element_text(face="bold", size=40),
                      axis.title.x = element_text(size=20, angle=0, face="bold"),
                      axis.title.y = element_text(size=25, face="bold"),
                      axis.text.y = element_text(size=13, face="bold"),
                      axis.text.x = element_text(size=13, face="bold"),
                      axis.ticks.y = element_blank(),
                      panel.grid.major = element_blank(),
                      panel.grid.minor = element_blank(),
                      panel.background = element_blank(),
                      legend.title=element_blank(),
                      legend.position = "none"
                )
hist <- axis_canvas(vmt.quad.plot, axis = "x") + 
        geom_histogram(data = vmt, aes(x = tmax), bins=30, fill='#FB8C00', color='gray60', alpha=0.5, size=0.1)
vmt.quad.plot <- insert_xaxis_grob(vmt.quad.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
vmt.quad.plot <- ggdraw(vmt.quad.plot)

## upt
upt[, tmax2 := tmax^2]
upt.quad <- felm(log_trips ~ tmax + tmax2
                  + trange + prcp.days + cloud + humid + wind 
                  | uza:month + uza:year 
                  | 0 
                  | uza
                  , data=upt)

t <- seq(from=-10, to=40, by=0.001)
d.t <- upt.quad$coef[1]*t + upt.quad$coef[2]*t^2
d.t <- d.t - max(d.t)
upt.t.data <- as.data.table(cbind(t, d.t))

upt.t.plot <- ggplot(upt.t.data) + 
    ylab("Log(PTT)") +
    xlab("Max. Temp. in \u2103") +
    geom_hline(yintercept = 0, linetype = 2, color="gray60") +
    geom_line(aes(x=t, y=d.t), color="#19A1FC", size=1.5) + 
    scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
    coord_cartesian(xlim=c(-10,40), ylim = c(-0.3, 0.01)) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=20, face="bold"),
          axis.title.y = element_text(size=25, face="bold"),
          axis.text.y = element_text(size=13, face="bold"),
          axis.text.x = element_text(size=13, face="bold"),
          axis.ticks.y = element_line(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_line(),
          legend.title=element_blank(),
          legend.position = "none"
    )
hist <- axis_canvas(upt.t.plot, axis = "x") + 
        geom_histogram(data = upt, aes(x = tmax), bins=30, fill='#19A1FC', color='gray60', alpha=0.5, size=0.1)
upt.t.plot <- insert_xaxis_grob(upt.t.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
upt.t.plot <- ggdraw(upt.t.plot)

##  margn effect temp
upt.quad.vcv <- upt.quad$clustervcv # get cluster robust variance covariance matrices

# get coefficients
upt.quad.tcoef <- upt.quad$coef[1]
upt.quad.t2coef <- upt.quad$coef[2]

# calculate marginal effects
marg <- data.table(tmax = seq(from=-10, to=40, by=0.001))
marg[, upt.quad.delta := upt.quad.tcoef + 2*upt.quad.t2coef*tmax]

# compute variance of marginal effect
marg[, upt.quad.se := sqrt(upt.quad.vcv["tmax","tmax"] + 4*(tmax^2)*upt.quad.vcv["tmax2","tmax2"] + 4*tmax*upt.quad.vcv[1,2])]

# 95% confidence bounds
marg[, upt.quad.upr := upt.quad.delta + 1.96*upt.quad.se]
marg[, upt.quad.lwr := upt.quad.delta - 1.96*upt.quad.se]

upt.quad.plot <- ggplot(data=marg) + 
                ylab("Marginal Effect of +1\u2103\non Log(PTT)") +
                xlab("Max. Temp. in \u2103") +
                geom_hline(aes(yintercept=0), linetype=2, color="gray60") +
                geom_ribbon(aes(x=tmax, ymin=upt.quad.lwr, ymax=upt.quad.upr), fill="gray60", alpha=0.25) + 
                geom_line(aes(x=tmax, y=upt.quad.delta), color="#19A1FC", size=1.5) +
                scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
                coord_cartesian(ylim=c(-0.003, 0.015)) +
                theme(plot.title = element_text(face="bold", size=40),
                      axis.title.x = element_text(size=20, angle=0, face="bold"),
                      axis.title.y = element_text(size=25, face="bold"),
                      axis.text.y = element_text(size=13, face="bold"),
                      axis.text.x = element_text(size=13, face="bold"),
                      axis.ticks.y = element_blank(),
                      panel.grid.major = element_blank(),
                      panel.grid.minor = element_blank(),
                      panel.background = element_blank(),
                      legend.title=element_blank(),
                      legend.position = "none"
                )
hist <- axis_canvas(upt.quad.plot, axis = "x") + 
        geom_histogram(data = upt, aes(x = tmax), bins=30, fill='#19A1FC', color='gray60', alpha=0.5, size=0.1)
upt.quad.plot <- insert_xaxis_grob(upt.quad.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
upt.quad.plot <- ggdraw(upt.quad.plot)

# plot----
png(filename = './robustness_checks.png', width=20, height=10, units = "in", type="cairo", res=300)
grid.arrange(vmt.main.bin.plot, 
             vmt.mon.bin.plot,
             vmt.t.plot,
             vmt.quad.plot,
             upt.main.bin.plot,
             upt.mon.bin.plot,
             upt.t.plot,
             upt.quad.plot,
             ncol=4, 
             nrow=2)
grid.text("a", x=unit(0.065, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("b", x=unit(0.315, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("c", x=unit(0.56, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("d", x=unit(0.84, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("e", x=unit(0.065, "npc"), y=unit(0.475, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("f", x=unit(0.31, "npc"), y=unit(0.475, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("g", x=unit(0.555, "npc"), y=unit(0.475, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("h", x=unit(0.83, "npc"), y=unit(0.475, "npc"), gp=gpar(fontsize=30,font=2))
dev.off()
```

```{r robust, out.width="100%", fig.cap = paste0(sfig_num("robust", display="cite"),": Plot of various robustness checks of main results. This figure plots the results of robustness checks that vary model specification for each of VMT and PTT. Panel (a) plots the results of pulling in the full distribution of temperature during each calendar month and estimating the regression using these five-degree Celsius binned counts of daily temperature. In this model the y-axis can be interpreted as the effect of one additional day in a specific temperature range out of the month on the log of VMT. Panel (b) also conducts a non-parametric binned regression, where average monthly temperatures are separated into five-degree Celsius bins. The y-axis in this regression can be interpreted as the effect of a monthly temperature falling into a particular temperature range. Panel (c) plots the fitted values from using a quadratic parameterization of monthly maximum temperatures. Panel (d) plots the marginal effects and confidence intervals of these marginal effects for the quadratic results presented in (c). Panels (e) through (h) replicate these modeling procedures for the PTT outcome variable. Across each of these parametric and non-parametric specifications of temperature, the main results from our splined regressions in the main text hold. Warmer temperatures, up to approximately 30C, increase both VMT and PTT. Error lines and regions represent 95\\% confidence intervals.")}
knitr::include_graphics("./robustness_checks.png")
```

```{r}
robust <- sfig_num(name = "robust", caption = "")
```

```{r vmt.monthlyproj, eval=FALSE}
# vehicle miles traveled by-state climate projection----
## vmt
# vmt_projection_spline <- function(data, dv, iv, lhs) {
#   ## keep spline knots with the fitted model results for plotting
#   spl_knots <- wtd.quantile(iv[!is.na(dv)], weights=dv[!is.na(dv)], probs = c(0.25, 0.5, 0.75), na.rm=T) # get quantiles of iv weighted by dv (cannot have missing values)
#   spl_boundary_knots <- range(iv, na.rm=T) # get range of independent variable
#   xs <- as.data.table(ns(iv, knots = spl_knots, Boundary.knots = spl_boundary_knots)) # get cubic splines of iv
#   setnames(xs, paste0("ns_", 1:ncol(xs))) # set names of splines
#   data <- cbind(data, xs)
#   ## set up felm model
#   rhs <- paste(c(names(xs), "trange", "prcp.days", "cloud", "humid", "wind"), collapse = " + ")
#   fe <- paste(c("state:month", "state:year"), collapse = " + ")
#   fmla <- as.formula(sprintf("%s ~ %s | %s | 0 | state", lhs, rhs, fe))
#   fit_felm <- felm(fmla, data=data)
#   result <- foreach::foreach (j = c('tmax_rcp26','tmax_rcp45','tmax_rcp85'), .combine = cbind) %do% {
#             pred_tmax <- state_nex[, get(j)] # prediction iv
#             xp <- as.data.table(ns(pred_tmax, knots = spl_knots, Boundary.knots = spl_boundary_knots))
#             result <- foreach::foreach (i = 1:nrow(xp), .combine = rbind) %do% {
#               this.x <- pred_tmax[i]
#               this.row <- c(unlist(xp[i, ]), 0, 0, 0, 0, 0) # include zeros for meteorological controls
#               temp.dt <- as.data.table(svycontrast(fit_felm, this.row))
#               temp.dt[, x := this.x]
#               setnames(temp.dt, c(paste0(j,"_coef"), paste0(j,"_se"), paste0(j,"_x")))
#               return(temp.dt)
#             }
#     return(result)
#   }
#   return(result)
# } # modified function from patrick baylis
# projection <- vmt_projection_spline(data=vmt, dv=vmt$log_vmt, iv=vmt$tmax, lhs="log_vmt")[, b:= 0] # main
# projection <- cbind(state_nex, projection)
# projection[, c('b', 'tmax_rcp85_x', 'tmax_rcp45_x', 'tmax_rcp26_x') := NULL] # remove unneeded vars
# projection[, ':=' (tmax_rcp26_coef = tmax_rcp26_coef - max(tmax_rcp26_coef),
#                    tmax_rcp45_coef = tmax_rcp45_coef - max(tmax_rcp45_coef),
#                    tmax_rcp85_coef = tmax_rcp85_coef - max(tmax_rcp85_coef))] # normalize so max==0
# projection[, tmax_rcp26_upr := tmax_rcp26_coef + 1.96*tmax_rcp26_se]
# projection[, tmax_rcp26_lwr := tmax_rcp26_coef - 1.96*tmax_rcp26_se]
# projection[, tmax_rcp26_se := NULL]
# projection[, tmax_rcp45_upr := tmax_rcp45_coef + 1.96*tmax_rcp45_se]
# projection[, tmax_rcp45_lwr := tmax_rcp45_coef - 1.96*tmax_rcp45_se]
# projection[, tmax_rcp45_se := NULL]
# projection[, tmax_rcp85_upr := tmax_rcp85_coef + 1.96*tmax_rcp85_se]
# projection[, tmax_rcp85_lwr := tmax_rcp85_coef - 1.96*tmax_rcp85_se]
# projection[, tmax_rcp85_se := NULL]
# saveRDS(projection, "./state_projection_bootstrap.rds")
# rm(projection)
state_projection <- readRDS("./state_projection_bootstrap.rds")

# calculate and plot by-month cumulative projection for rcp26/45/85----
## need to use 2006-2020 tmax baseline and calculate differences by state then cumulate
mean2006_2020 <- state_projection[year(monthyr) >= 2006 & year(monthyr) < 2020, 
                                  .(base_rcp26_coef = mean(tmax_rcp26_coef),
                                    base_rcp45_coef = mean(tmax_rcp45_coef),
                                    base_rcp85_coef = mean(tmax_rcp85_coef)), 
                                  by=.(month(monthyr), state)]
state_projection[, month := month(monthyr)]
setkey(state_projection, state, month)
setkey(mean2006_2020, state, month)
state_projection <- merge(state_projection, mean2006_2020)
rm(mean2006_2020)

## get baseline vmt (millions) by state for 2013-2017 (most recent full five-year period)
baseline_vmt <- vmt[year(monthyr) >= 2013 & year(monthyr) <= 2017, .(vmt = mean(vmt, na.rm=T)), by=.(state, month(monthyr))]
setkey(baseline_vmt, state, month)
setkey(state_projection, state, month)
state_projection <- merge(state_projection, baseline_vmt)
rm(baseline_vmt)

## calculate changes relative to baseline
## exponentiate coefficients to get percentage change difference from fixed baseline
state_projection[, ':=' (delta_rcp26 = 100*(exp(tmax_rcp26_coef) - exp(base_rcp26_coef)),
                         delta_rcp45 = 100*(exp(tmax_rcp45_coef) - exp(base_rcp45_coef)),
                         delta_rcp85 = 100*(exp(tmax_rcp85_coef) - exp(base_rcp85_coef)))] 

## keep only year 2020+
state_projection <- state_projection[year(monthyr) >= 2020]

## get by-month
month20402050 <- state_projection[year(monthyr) >= 2040 & year(monthyr) <= 2050, 
                                 .(delta_rcp26 = mean(delta_rcp26, na.rm=T),
                                   delta_rcp45 = mean(delta_rcp45, na.rm=T),
                                   delta_rcp85 = mean(delta_rcp85, na.rm=T)),
                                 by=.(month, state=as.character(state))
                                 ]
month20402050 <- melt(month20402050, id.vars = c('month', 'state'), variable.name = "rcp", value.name = 'vmt')

month20902100 <- state_projection[year(monthyr) >= 2090 & year(monthyr) <= 2100, 
                                 .(delta_rcp26 = mean(delta_rcp26, na.rm=T),
                                   delta_rcp45 = mean(delta_rcp45, na.rm=T),
                                   delta_rcp85 = mean(delta_rcp85, na.rm=T)),
                                 by=.(month, state=as.character(state))
                                 ]
month20902100 <- melt(month20902100, id.vars = c('month','state'), variable.name = "rcp", value.name = 'vmt')

## plot data
names(states)[6] <- "state" # fix names
states$state <- as.character(states$state)
month2050rcp26 <- month20402050[rcp=="delta_rcp26", .(vmt), by=.(month, state)]
month2099rcp26 <- month20902100[rcp=="delta_rcp26", .(vmt), by=.(month, state)]
month2050rcp26 <- dcast(month2050rcp26, state ~ paste0('month',month), value.var="vmt")
month2099rcp26 <- dcast(month2099rcp26, state ~ paste0('month',month),, value.var="vmt")
month2050rcp26 <- sp::merge(states, month2050rcp26, by="state", duplicateGeoms=TRUE)
month2099rcp26 <- sp::merge(states, month2099rcp26, by="state", duplicateGeoms=TRUE)

month2050rcp45 <- month20402050[rcp=="delta_rcp45", .(vmt), by=.(month, state)]
month2099rcp45 <- month20902100[rcp=="delta_rcp45", .(vmt), by=.(month, state)]
month2050rcp45 <- dcast(month2050rcp45, state ~ paste0('month',month),, value.var="vmt")
month2099rcp45 <- dcast(month2099rcp45, state ~ paste0('month',month),, value.var="vmt")
month2050rcp45 <- sp::merge(states, month2050rcp45, by="state", duplicateGeoms=TRUE)
month2099rcp45 <- sp::merge(states, month2099rcp45, by="state", duplicateGeoms=TRUE)

month2050rcp85 <- month20402050[rcp=="delta_rcp85", .(vmt), by=.(month, state)]
month2099rcp85 <- month20902100[rcp=="delta_rcp85", .(vmt), by=.(month, state)]
month2050rcp85 <- dcast(month2050rcp85, state ~ paste0('month',month),, value.var="vmt")
month2099rcp85 <- dcast(month2099rcp85, state ~ paste0('month',month),, value.var="vmt")
month2050rcp85 <- sp::merge(states, month2050rcp85, by="state", duplicateGeoms=TRUE)
month2099rcp85 <- sp::merge(states, month2099rcp85, by="state", duplicateGeoms=TRUE)

# rcp26----
p1 <- spplot(
       obj=month2050rcp26,
       zcol=c('month1','month2','month3','month4','month5','month6','month7','month8','month9','month10','month11','month12'),         
       names.attr = paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December")), # set names for strip text
       strip=TRUE, # include strip text
       col.regions=colorRampPalette(c("#bcd168", "#d5ed76", "#ffc67f", "#FB8C00", "#e27e00", "#9e3f00"))(100), # set colors
       as.table = TRUE,
       # col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2), height=0.5), # make legend on bottom of plot
       at=seq(min(c(month2099rcp85$month6,month2099rcp85$month7,month2099rcp85$month8,month2099rcp85$month9), na.rm=TRUE), max(c(month2099rcp85$month11,month2099rcp85$month12,month2099rcp85$month1,month2099rcp85$month2), na.rm=TRUE)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       layout=c(4,3) # horizontal layout
       )
p1 <- update(p1, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p1 <- p1 + layer_(sp.lines(obj=us, fill="white", col="black"))
p1 <- p1 + layer(sp.lines(obj=states,
                              fill="transparent",
                              col="black",
                              cex=0.1))

p2 <- spplot(
       obj=month2099rcp26,
       zcol=c('month1','month2','month3','month4','month5','month6','month7','month8','month9','month10','month11','month12'),         
       names.attr = paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December")), # set names for strip text
       strip=TRUE, # include strip text
       col.regions=colorRampPalette(c("#bcd168", "#d5ed76", "#ffc67f", "#FB8C00", "#e27e00", "#9e3f00"))(100), # set colors
       as.table = TRUE,
       # col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2), height=0.5), # make legend on bottom of plot
       at=seq(min(c(month2099rcp85$month6,month2099rcp85$month7,month2099rcp85$month8,month2099rcp85$month9), na.rm=TRUE), max(c(month2099rcp85$month11,month2099rcp85$month12,month2099rcp85$month1,month2099rcp85$month2), na.rm=TRUE)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       layout=c(4,3) # horizontal layout
       )
p2 <- update(p2, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p2 <- p2 + layer_(sp.lines(obj=us, fill="white", col="black"))
p2 <- p2 + layer(sp.lines(obj=states,
                              fill="transparent",
                              col="black",
                              cex=0.1))

## get legend for plotting
key <- draw.colorkey(p1$legend[[1]]$args$key)
p1$legend <- NULL
p2$legend <- NULL
p1$sub <- NULL
p2$sub <- NULL

png(filename = './bymonth_rcp26.png', width=9.75, height=11, units = "in", type="cairo", res=300)
grid.arrange(rectGrob(gp=gpar(col="white")), 
             p1, 
             p2, 
             key, 
             textGrob("Percentage Change in VMT, RCP2.6", gp=gpar(fontsize=23, font=2)),
             nrow=5, 
             heights=c(0.02, 0.43, 0.43, 0.07, 0.05))
grid.text("2040-2050", x=unit(0.5, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=21, font=2))
grid.text("2090-2100", x=unit(0.5, "npc"), y=unit(0.545, "npc"), gp=gpar(fontsize=21, font=2))
dev.off()

# rcp45----
p1 <- spplot(
       obj=month2050rcp45,
       zcol=c('month1','month2','month3','month4','month5','month6','month7','month8','month9','month10','month11','month12'),         
       names.attr = paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December")), # set names for strip text
       strip=TRUE, # include strip text
       col.regions=colorRampPalette(c("#bcd168", "#d5ed76", "#ffc67f", "#FB8C00", "#e27e00", "#9e3f00"))(100), # set colors
       as.table = TRUE,
       # col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2), height=0.5), # make legend on bottom of plot
       at=seq(min(c(month2099rcp85$month6,month2099rcp85$month7,month2099rcp85$month8,month2099rcp85$month9), na.rm=TRUE), max(c(month2099rcp85$month11,month2099rcp85$month12,month2099rcp85$month1,month2099rcp85$month2), na.rm=TRUE)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       layout=c(4,3) # horizontal layout
       )
p1 <- update(p1, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p1 <- p1 + layer_(sp.lines(obj=us, fill="white", col="black"))
p1 <- p1 + layer(sp.lines(obj=states,
                              fill="transparent",
                              col="black",
                              cex=0.1))

p2 <- spplot(
       obj=month2099rcp45,
       zcol=c('month1','month2','month3','month4','month5','month6','month7','month8','month9','month10','month11','month12'),         
       names.attr = paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December")), # set names for strip text
       strip=TRUE, # include strip text
       col.regions=colorRampPalette(c("#bcd168", "#d5ed76", "#ffc67f", "#FB8C00", "#e27e00", "#9e3f00"))(100), # set colors
       as.table = TRUE,
       # col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2), height=0.5), # make legend on bottom of plot
       at=seq(min(c(month2099rcp85$month6,month2099rcp85$month7,month2099rcp85$month8,month2099rcp85$month9), na.rm=TRUE), max(c(month2099rcp85$month11,month2099rcp85$month12,month2099rcp85$month1,month2099rcp85$month2), na.rm=TRUE)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       layout=c(4,3) # horizontal layout
       )
p2 <- update(p2, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p2 <- p2 + layer_(sp.lines(obj=us, fill="white", col="black"))
p2 <- p2 + layer(sp.lines(obj=states,
                              fill="transparent",
                              col="black",
                              cex=0.1))

## get legend for plotting
key <- draw.colorkey(p1$legend[[1]]$args$key)
p1$legend <- NULL
p2$legend <- NULL
p1$sub <- NULL
p2$sub <- NULL

png(filename = './bymonth_rcp45.png', width=9.75, height=11, units = "in", type="cairo", res=300)
grid.arrange(rectGrob(gp=gpar(col="white")), 
             p1, 
             p2, 
             key, 
             textGrob("Percentage Change in VMT, RCP4.5", gp=gpar(fontsize=23, font=2)),
             nrow=5, 
             heights=c(0.02, 0.43, 0.43, 0.07, 0.05))
grid.text("2040-2050", x=unit(0.5, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=21, font=2))
grid.text("2090-2100", x=unit(0.5, "npc"), y=unit(0.545, "npc"), gp=gpar(fontsize=21, font=2))
dev.off()
```

```{r monthlyproj26, out.width="100%", fig.cap = paste0(sfig_num("monthlyproj26", display="cite"),": Plot of monthly projected mean percentage changes in vehicle miles traveled 2040-2050 and 2090-2100 across the RCP2.6 emissions scenario. This figure replicates the methods of Figure 4 from the main text across the RCP2.6 emissions scenario. The RCP2.6 emissions scenario results in lower projected changes in future VMT.")}
knitr::include_graphics("./bymonth_rcp26.png")
```

```{r}
monthlyproj26 <- sfig_num(name = "monthlyproj26", caption = "")
```

```{r monthlyproj45, out.width="100%", fig.cap = paste0(sfig_num("monthlyproj45", display="cite"),": Plot of monthly projected mean percentage changes in vehicle miles traveled 2040-2050 and 2090-2100 across the RCP4.5 emissions scenario. This figure replicates the methods of Figure 4 from the main text across the RCP4.5 emissions scenario. Winter months observe increases in VMT while summer months observe decreases. The RCP4.5 emissions scenario results in middle-range projected changes in future VMT.")}
knitr::include_graphics("./bymonth_rcp45.png")
```

```{r}
monthlyproj45 <- sfig_num(name = "monthlyproj45", caption = "")
```

```{r ptt.fig3, eval=FALSE}
load('./replication.RData')
# plot historical tmax histogram with indication of peak ptt----
dum <- data.frame(point = sort(rep(letters[1], 10)), x = runif(10), y = runif(10))
tmax.hist <- ggplot(data = upt, aes(x = tmax)) +
                  geom_histogram(aes(y=..density..), fill='#d5ed76', color='gray40', alpha=0.8, bins=30) +
                  ylab("Density") + 
                  xlab("Max. Temp. in \u2103") +
                  geom_vline(xintercept = upt.main[b==0 & coef_norm==max(coef_norm), x], color='#19A1FC', linetype=1, size=3) +
                  geom_line(data = dum, aes(color = point, x = x, y = y), alpha = 0, size=4) + 
                  scale_color_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                     limits=c("#19A1FC"), 
                     values=c("#19A1FC"),
                     labels=c("Peak\nHistorical\nPTT")) +
                  coord_cartesian(ylim=c(0, 0.045), xlim=c(-12,42)) +
                  theme(plot.title = element_text(face="bold", size=40),
                        axis.title.x = element_text(size=20, face="bold"),
                        axis.title.y = element_text(size=20, face="bold"),
                        axis.text.y = element_text(size=12, face="bold"),
                        axis.text.x = element_text(size=12, face="bold"),
                        panel.grid.major = element_blank(),
                        panel.grid.minor = element_blank(),
                        panel.background = element_blank(),
                        axis.line.x = element_line(),
                        axis.line.y = element_line(),
                        axis.ticks.y = element_line(),
                        legend.title=element_blank(),
                        legend.position = c(.01, 0.7),
                        legend.text=element_text(size=12, face="bold"),
                        legend.key = element_rect(colour = "black", fill="white")
                  )

# plot time series of monthly gcm projections by uza----
## calculate temperature anomalies relative to baseline
gcm_baseline <- uza_nex[year(monthyr) >= 2006 & year(monthyr) < 2020, .(baseline_rcp26 = mean(tmax_rcp26),
                        baseline_rcp45 = mean(tmax_rcp45),
                        baseline_rcp85 = mean(tmax_rcp85)),
                    by=.(month(monthyr), uza)]
uza_nex[, month := month(monthyr)]
setkey(gcm_baseline, uza, month)
setkey(uza_nex, uza, month)
gcm <- merge(uza_nex, gcm_baseline)
gcm <- gcm[year(monthyr) >= 2020,] # keep only 2020+
gcm[, ':=' (rcp26_anom = tmax_rcp26 - baseline_rcp26,
            rcp45_anom = tmax_rcp45 - baseline_rcp45,
            rcp85_anom = tmax_rcp85 - baseline_rcp85)]

dum <- data.frame(point = sort(rep(letters[1:3], 10)), x = gcm$monthyr, y = runif(30))
gcm_plot <- ggplot(gcm, aes(x=monthyr, group=uza)) + 
                  geom_hline(aes(yintercept=0), linetype=2, color='gray80') +
                  geom_line(aes(y=rcp26_anom), color='#d5ed76', size=0.02, alpha=0.3) + 
                  geom_line(aes(y=rcp45_anom), color='gray60', size=0.02, alpha=0.3) + 
                  geom_line(aes(y=rcp85_anom), color='#19A1FC', size=0.02, alpha=0.3) +
                  ylab("Max. Temp. Anomaly in \u2103") +
                  xlab("Calendar Month") +
                  geom_line(data = dum, aes(color = point, x = x, y = y, group=NULL), alpha = 0, size=4) + 
                  scale_color_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                     limits=c("#d5ed76","gray60","#19A1FC"), 
                     values=c("#d5ed76","gray60","#19A1FC"),
                     labels=c("RCP2.6","RCP4.5","RCP8.5")) +
                  coord_cartesian(xlim=c(2020,2100)) +
                  theme(plot.title = element_text(face="bold", size=40),
                        axis.title.x = element_text(size=20, face="bold"),
                        axis.title.y = element_text(size=20, face="bold"),
                        axis.text.y = element_text(size=12, face="bold"),
                        axis.text.x = element_text(size=12, face="bold"),
                        panel.grid.major = element_blank(),
                        panel.grid.minor = element_blank(),
                        panel.background = element_blank(),
                        axis.line.x = element_line(),
                        axis.line.y = element_line(),
                        axis.ticks.y = element_line(),
                        legend.title=element_blank(),
                        legend.position = c(.01, 0.7),
                        legend.text=element_text(size=12, face="bold"),
                        legend.key = element_rect(colour = "black", fill="white")
                  )

# public transit trips by-uza climate projection----
## ptt
# ptt_projection_spline <- function(data, dv, iv, lhs) {
#   ## keep spline knots with the fitted model results for plotting
#   spl_knots <- wtd.quantile(iv[!is.na(dv)], weights=dv[!is.na(dv)], probs = c(0.25, 0.5, 0.75), na.rm=T) # get quantiles of iv weighted by dv (cannot have missing values)
#   spl_boundary_knots <- range(iv, na.rm=T) # get range of independent variable
#   xs <- as.data.table(ns(iv, knots = spl_knots, Boundary.knots = spl_boundary_knots)) # get cubic splines of iv
#   setnames(xs, paste0("ns_", 1:ncol(xs))) # set names of splines
#   data <- cbind(data, xs)
#   ## set up felm model
#   rhs <- paste(c(names(xs), "trange", "prcp.days", "cloud", "humid", "wind"), collapse = " + ")
#   fe <- paste(c("uza:month", "uza:year"), collapse = " + ")
#   fmla <- as.formula(sprintf("%s ~ %s | %s | 0 | uza", lhs, rhs, fe))
#   fit_felm <- felm(fmla, data=data)
#   result <- foreach::foreach (j = c('tmax_rcp26','tmax_rcp45','tmax_rcp85'), .combine = cbind) %dopar% {
#             pred_tmax <- uza_nex[, get(j)] # prediction iv
#             xp <- as.data.table(ns(pred_tmax, knots = spl_knots, Boundary.knots = spl_boundary_knots))
#             result <- foreach::foreach (i = 1:nrow(xp), .combine = rbind) %do% {
#               this.x <- pred_tmax[i]
#               this.row <- c(unlist(xp[i, ]), 0, 0, 0, 0, 0) # include zeros for meteorological controls
#               temp.dt <- as.data.table(svycontrast(fit_felm, this.row))
#               temp.dt[, x := this.x]
#               setnames(temp.dt, c(paste0(j,"_coef"), paste0(j,"_se"), paste0(j,"_x")))
#               return(temp.dt)
#             }
#     return(result)
#   }
#   return(result)
# } # modified function from patrick baylis
# projection <- ptt_projection_spline(data=upt, dv=upt$log_trips, iv=upt$tmax, lhs="log_trips")[, b:= 0] # main
# projection <- cbind(uza_nex, projection)
# projection[, c('b', 'tmax_rcp85_x', 'tmax_rcp45_x', 'tmax_rcp26_x') := NULL] # remove unneeded vars
# projection[, ':=' (tmax_rcp26_coef = tmax_rcp26_coef - max(tmax_rcp26_coef),
#                    tmax_rcp45_coef = tmax_rcp45_coef - max(tmax_rcp45_coef),
#                    tmax_rcp85_coef = tmax_rcp85_coef - max(tmax_rcp85_coef))] # normalize so max==0
# projection[, tmax_rcp26_upr := tmax_rcp26_coef + 1.96*tmax_rcp26_se]
# projection[, tmax_rcp26_lwr := tmax_rcp26_coef - 1.96*tmax_rcp26_se]
# projection[, tmax_rcp26_se := NULL]
# projection[, tmax_rcp45_upr := tmax_rcp45_coef + 1.96*tmax_rcp45_se]
# projection[, tmax_rcp45_lwr := tmax_rcp45_coef - 1.96*tmax_rcp45_se]
# projection[, tmax_rcp45_se := NULL]
# projection[, tmax_rcp85_upr := tmax_rcp85_coef + 1.96*tmax_rcp85_se]
# projection[, tmax_rcp85_lwr := tmax_rcp85_coef - 1.96*tmax_rcp85_se]
# projection[, tmax_rcp85_se := NULL]
# saveRDS(projection, "./uza_projection_bootstrap.rds")
# rm(projection)
uza_projection <- readRDS("./uza_projection_bootstrap.rds")

# calculate and plot by-month cumulative projection for rcp26/45/85----
## need to use 2006-2020 tmax and calculate differences by state then cumulate
mean2006_2020 <- uza_projection[year(monthyr) >= 2006 & year(monthyr) < 2020, 
                                  .(base_rcp26_coef = mean(tmax_rcp26_coef),
                                    base_rcp45_coef = mean(tmax_rcp45_coef),
                                    base_rcp85_coef = mean(tmax_rcp85_coef)), 
                                  by=.(month(monthyr), uza)]
uza_projection[, month := month(monthyr)]
setkey(uza_projection, uza, month)
setkey(mean2006_2020, uza, month)
uza_projection <- merge(uza_projection, mean2006_2020)
rm(mean2006_2020)

## get baseline ptt by uza for 2008-2017 (most recent ten-year period)
## use ten year period to keep greater number of cities
baseline_ptt <- upt[year(monthyr) >= 2008 & year(monthyr) <= 2017, .(trips = mean(trips, na.rm=T)), by=.(uza = as.integer(as.character(uza)), month(monthyr))]
setkey(baseline_ptt, uza, month)
setkey(uza_projection, uza, month)
uza_projection <- merge(uza_projection, baseline_ptt)
rm(baseline_ptt)

## calculate changes relative to baseline
## exponentiate coefficients to get percentage change, convert to trips, difference from fixed baseline
uza_projection[, ':=' (delta_rcp26 = exp(tmax_rcp26_coef)*trips - exp(base_rcp26_coef)*trips,
                         delta_rcp45 = exp(tmax_rcp45_coef)*trips - exp(base_rcp45_coef)*trips,
                         delta_rcp85 = exp(tmax_rcp85_coef)*trips - exp(base_rcp85_coef)*trips)] 

## keep only year 2020+
uza_projection <- uza_projection[year(monthyr) >= 2020]

## cumulate changes by uza across month-years
setkey(uza_projection, uza, monthyr)
uza_cumul_changes <- uza_projection[, .(cumul_rcp26 = cumsum(delta_rcp26),
                                            cumul_rcp45 = cumsum(delta_rcp45),
                                            cumul_rcp85 = cumsum(delta_rcp85),
                                          monthyr = unique(monthyr)), by=.(uza)]

## cumulate changes by full us
us_cumul_changes <- uza_cumul_changes[, .(cumul_rcp26 = sum(cumul_rcp26),
                                            cumul_rcp45 = sum(cumul_rcp45),
                                            cumul_rcp85 = sum(cumul_rcp85)), 
                                        by=.(monthyr)]

dum <- data.frame(point = sort(rep(letters[1:3], 10)), x = us_cumul_changes$monthyr, y = runif(30))
us_cumul_plot <- ggplot(us_cumul_changes, aes(x=monthyr)) + 
                  geom_hline(aes(yintercept=0), linetype=2, color='gray80') +
                  geom_line(aes(y=cumul_rcp26/1000000000), color='#d5ed76', size=3) + 
                  geom_line(aes(y=cumul_rcp45/1000000000), color='gray60', size=3) + 
                  geom_line(aes(y=cumul_rcp85/1000000000), color='#19A1FC', size=3) +
                  ylab("Cumulative Added PTT in Billions of Trips") +
                  xlab("Calendar Month") +
                  geom_line(data = dum, aes(color = point, x = x, y = y), alpha = 0, size=4) + 
                  scale_color_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                     limits=c("#d5ed76","gray60","#19A1FC"), 
                     values=c("#d5ed76","gray60","#19A1FC"),
                     labels=c("RCP2.6","RCP4.5","RCP8.5")) +
                  coord_cartesian(ylim=c(0,6)) +
                  # scale_y_continuous(breaks=c(0, 250, 500, 750, 1000)) +
                  theme(plot.title = element_text(face="bold", size=40),
                        axis.title.x = element_text(size=20, face="bold"),
                        axis.title.y = element_text(size=25, face="bold"),
                        axis.text.y = element_text(size=12, face="bold"),
                        axis.text.x = element_text(size=12, face="bold"),
                        panel.grid.major = element_blank(),
                        panel.grid.minor = element_blank(),
                        panel.background = element_blank(),
                        axis.line.x = element_line(),
                        axis.line.y = element_line(),
                        axis.ticks.y = element_line(),
                        legend.title=element_blank(),
                        legend.position = c(.2, 0.7),
                        legend.text=element_text(size=17, face="bold"),
                        legend.key = element_rect(colour = "black", fill="white")
                  )

# calculate and plot 2040-2050 and 2090-2100 by-month projections----
month20402050 <- uza_projection[year(monthyr) >= 2040 & year(monthyr) <= 2050, 
                                 .(delta_rcp26 = sum(delta_rcp26),
                                   delta_rcp45 = sum(delta_rcp45),
                                   delta_rcp85 = sum(delta_rcp85)),
                                 by=.(month)
                                 ]
month20402050 <- melt(month20402050, id.vars = 'month', variable.name = "rcp", value.name = 'trips')

month20902100 <- uza_projection[year(monthyr) >= 2090 & year(monthyr) <= 2100, 
                                 .(delta_rcp26 = sum(delta_rcp26),
                                   delta_rcp45 = sum(delta_rcp45),
                                   delta_rcp85 = sum(delta_rcp85)),
                                 by=.(month)
                                 ]
month20902100 <- melt(month20902100, id.vars = 'month', variable.name = "rcp", value.name = 'trips')

## plot
rcp.colors <- c(delta_rcp26 = "#d5ed76", delta_rcp45 = "gray60", delta_rcp85 ="#19A1FC") # set colors

month20402050.plot <- ggplot(data=month20402050) +
    geom_hline(yintercept = 0, color="gray60") +
    geom_col(aes(x=month, y=trips/1000000, width=0.8, fill=factor(rcp)), color="gray60", position = "dodge") +
    coord_cartesian(ylim = c(-20,200), xlim = c(0.5,12.5)) +
    scale_x_continuous(breaks=c(1:12), labels=c('J','F','M','A','M','J','J','A','S','O','N','D')) +
    ylab("Added PTT in Millions") +
    xlab("Months, 2040-2050") +
    scale_fill_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                     values=rcp.colors,
                     labels=c("RCP2.6","RCP4.5","RCP8.5")) +
    theme(plot.title = element_text(face="bold", size=40),
                        axis.title.x = element_text(size=20, face="bold"),
                        axis.title.y = element_blank(),
                        axis.text.y = element_blank(),
                        axis.text.x = element_text(size=12, face="bold"),
                        panel.grid.major = element_blank(),
                        panel.grid.minor = element_blank(),
                        panel.background = element_blank(),
                        axis.line.x = element_line(),
                        axis.line.y = element_blank(),
                        axis.ticks.y = element_blank(),
                        legend.title=element_blank(),
                        legend.position = c(.7, 0.6),
                        legend.text=element_text(size=12, face="bold"),
                        legend.key = element_rect(colour = "black", fill="white")
                  )

month209002100.plot <- ggplot(data=month20902100) +
    geom_hline(yintercept = 0, color="gray60") +
    geom_col(aes(x=month, y=trips/1000000, width=0.8, fill=factor(rcp)), color="gray60", position = "dodge") +
    coord_cartesian(ylim = c(-20,200), xlim = c(0.5,12.5)) +
    scale_x_continuous(breaks=c(1:12), labels=c('J','F','M','A','M','J','J','A','S','O','N','D')) +
    ylab("Added PTT in Millions") +
    xlab("Months, 2090-2100") +
    scale_fill_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                     values=rcp.colors,
                     labels=c("RCP2.6","RCP4.5","RCP8.5")) +
    theme(plot.title = element_text(face="bold", size=40),
                        axis.title.x = element_text(size=20, face="bold"),
                        axis.title.y = element_blank(),
                        axis.text.y = element_blank(),
                        axis.text.x = element_text(size=12, face="bold"),
                        panel.grid.major = element_blank(),
                        panel.grid.minor = element_blank(),
                        panel.background = element_blank(),
                        axis.line.x = element_line(),
                        axis.line.y = element_blank(),
                        axis.ticks.y = element_blank(),
                        legend.title=element_blank(),
                        legend.position = "none"
                  )

justy.month <- ggplot(data=month20902100) +
    geom_hline(yintercept = 0, color="gray60") +
    geom_col(aes(x=month, y=trips/1000000, width=0.8, fill=factor(rcp)), color="gray60", position = "dodge") +
    coord_cartesian(ylim = c(-20,200), xlim = c(0.5,12.5)) +
    scale_x_continuous(breaks=c(1:12), labels=c('J','F','M','A','M','J','J','A','S','O','N','D')) +
    ylab("Added PTT in Millions") +
    xlab("Months, 2090-2100") +
    scale_fill_manual(guide = guide_legend(override.aes = list(alpha = 1)),
                     values=rcp.colors,
                     labels=c("RCP2.6","RCP4.5","RCP8.5")) +
    theme(plot.title = element_text(face="bold", size=40),
                        axis.title.x = element_text(size=20, face="bold"),
                        axis.title.y = element_text(size=20, face="bold"),
                        axis.text.y = element_text(size=12, face="bold"),
                        axis.text.x = element_text(size=12, face="bold"),
                        panel.grid.major = element_blank(),
                        panel.grid.minor = element_blank(),
                        panel.background = element_blank(),
                        axis.line.x = element_line(),
                        axis.line.y = element_line(),
                        axis.ticks.y = element_line(),
                        legend.title=element_blank(),
                        legend.position = c(.05, 0.2),
                        legend.text=element_text(size=16, face="bold"),
                        legend.key = element_rect(colour = "black", fill="white")
                  )
justy.month <- gtable_filter(ggplotGrob(justy.month), 'axis-l|ylab', trim=F)
justy.month <- ggdraw(justy.month)

# calculate and plot 2040-2050 and 2090-2100 by-city cumulative map projections----
uza_projection <- readRDS("./uza_projection_bootstrap.rds")

## need to merge in state ids based on uza centroids
names(states)[6] <- "state" # fix names
states$state <- as.character(states$state)
points <- uza_projection[, .(lat = unique(lat), lon = unique(lon)), by=.(uza)]
coordinates(points) <- ~lon + lat
projection(points) <- projection(states)
uza.state <- cbind(as.data.table(points), as.data.table(over(points, states)))
uza.state <- uza.state[, .(uza, state)]
setkey(uza.state, uza)
setkey(uza_projection, uza)
uza_projection <- merge(uza_projection, uza.state)
uza_projection[uza==13, state := "California"]
uza_projection[uza==17, state := "Florida"]
uza_projection[uza==34, state := "Virginia"]

## average uza_projection to state-level
uza_projection <- uza_projection[, lapply(.SD, mean), .SDcols=c("tmax_rcp45","tmax_rcp85","tmax_rcp26","tmax_rcp26_coef","tmax_rcp45_coef","tmax_rcp85_coef","tmax_rcp26_upr","tmax_rcp26_lwr","tmax_rcp45_upr","tmax_rcp45_lwr","tmax_rcp85_upr","tmax_rcp85_lwr"),
                                 by=.(state, monthyr, month)]

## need to use 2006-2020 tmax baseline and calculate differences by state then cumulate
mean2006_2020 <- uza_projection[year(monthyr) >= 2006 & year(monthyr) < 2020, 
                                  .(base_rcp26_coef = mean(tmax_rcp26_coef),
                                    base_rcp45_coef = mean(tmax_rcp45_coef),
                                    base_rcp85_coef = mean(tmax_rcp85_coef)), 
                                  by=.(month(monthyr), state)]
uza_projection[, month := month(monthyr)]
setkey(uza_projection, state, month)
setkey(mean2006_2020, state, month)
uza_projection <- merge(uza_projection, mean2006_2020)
rm(mean2006_2020)

## get baseline ptt by state for 2008-2017 (most recent full ten-year period)
setkey(uza.state, uza)
setkey(upt, uza)
upt <- merge(upt, uza.state, all.x=T)
baseline_ptt <- upt[, .(trips = sum(trips, na.rm=T)), 
                    by=.(state, monthyr)] # sum to state-level
baseline_ptt <- baseline_ptt[year(monthyr) >= 2008 & year(monthyr) <= 2017, 
                            .(trips = mean(trips, na.rm=T)), 
                            by=.(state, month = month(monthyr))]
setkey(baseline_ptt, state, month)
setkey(uza_projection, state, month)
uza_projection <- merge(uza_projection, baseline_ptt)
rm(baseline_ptt)

## calculate changes relative to baseline
## exponentiate coefficients to get percentage change, convert to trips, difference from fixed baseline
uza_projection[, ':=' (delta_rcp26 = exp(tmax_rcp26_coef)*trips - exp(base_rcp26_coef)*trips,
                         delta_rcp45 = exp(tmax_rcp45_coef)*trips - exp(base_rcp45_coef)*trips,
                         delta_rcp85 = exp(tmax_rcp85_coef)*trips - exp(base_rcp85_coef)*trips)] 

## keep only year 2020+
uza_projection <- uza_projection[year(monthyr) >= 2020]

## cumulate changes by state across month-years
setkey(uza_projection, state, monthyr)
uza_cumul_changes <- uza_projection[, .(cumul_rcp26 = cumsum(delta_rcp26),
                                            cumul_rcp45 = cumsum(delta_rcp45),
                                            cumul_rcp85 = cumsum(delta_rcp85),
                                          monthyr = unique(monthyr)), by=.(state)]

## mean of cumulative changes by state in 2050, 2099 by model
state20502099 <- uza_cumul_changes[year(monthyr)==2050 | year(monthyr)==2099, 
                                     .(cumul_rcp26 = mean(cumul_rcp26),
                                       cumul_rcp45 = mean(cumul_rcp45),
                                       cumul_rcp85 = mean(cumul_rcp85)),
                                     by=.(year(monthyr), state)]

## plot
## ptt
names(states)[6] <- "state" # fix names
states$state <- as.character(states$state)
cumul.ptt.2050 <- state20502099[year==2050, 
                                .(cumul_rcp26 = log(cumul_rcp26),
                                 cumul_rcp45 = log(cumul_rcp45),
                                 cumul_rcp85 = log(cumul_rcp85)
                                     ), by=.(state)]
cumul.ptt.2050[, state := as.character(state)]
cumul.ptt.2050 <- sp::merge(states, cumul.ptt.2050, by="state")

cumul.ptt.2099 <- state20502099[year==2099, 
                                .(cumul_rcp26 = log(cumul_rcp26),
                                 cumul_rcp45 = log(cumul_rcp45),
                                 cumul_rcp85 = log(cumul_rcp85)
                                     ), by=.(state)]
cumul.ptt.2099[, state := as.character(state)]
cumul.ptt.2099 <- sp::merge(states, cumul.ptt.2099, by="state")

p1 <- spplot(
       obj=cumul.ptt.2050,
       zcol=c("cumul_rcp26","cumul_rcp45","cumul_rcp85"),
       strip=FALSE,
       col.regions=colorRampPalette(c("#ddf1ff", "#c4e7ff", "#87cdff","#19A1FC", "#0074c1"))(100), # set colors
       as.table = TRUE,
       col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2), height=0.5), # make legend on bottom of plot
       at=seq(min(cumul.ptt.2050$cumul_rcp26, na.rm=TRUE), max(cumul.ptt.2099$cumul_rcp85, na.rm=TRUE)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       sub=list("log(Cumulative Added Public Transit Trips)", cex=1.4),
       layout=c(3,1) # horizontal layout
       )
p1 <- update(p1, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p1 <- p1 + layer_(sp.lines(obj=us, fill="white", col="black"))
p1 <- p1 + layer(sp.lines(obj=states,
                              fill="transparent",
                              col="black",
                              cex=0.1))

p2 <- spplot(
       obj=cumul.ptt.2099,
       zcol=c("cumul_rcp26","cumul_rcp45","cumul_rcp85"),
       strip=FALSE,
       col.regions=colorRampPalette(c("#ddf1ff", "#c4e7ff", "#87cdff","#19A1FC", "#0074c1"))(100), # set colors
       as.table = TRUE,
       col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2), height=0.5), # make legend on bottom of plot
       at=seq(min(cumul.ptt.2050$cumul_rcp26, na.rm=TRUE), max(cumul.ptt.2099$cumul_rcp85, na.rm=TRUE)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       sub=list("log(Cumulative Added Public Transit Trips)", cex=1.4),
       layout=c(3,1) # horizontal layout
       )
p2 <- update(p2, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p2 <- p2 + layer_(sp.lines(obj=us, fill="white", col="black"))
p2 <- p2 + layer(sp.lines(obj=states,
                              fill="transparent",
                              col="black",
                              cex=0.1))

## get legend for plotting
key <- draw.colorkey(p1$legend[[1]]$args$key)
p1$legend <- NULL
p2$legend <- NULL
p1$sub <- NULL
p2$sub <- NULL

# plot----
png(filename = './ptt_projection.png', width=14, height=15, units = "in", type="cairo", res=300)
grid.arrange(arrangeGrob(rectGrob(gp=gpar(col="white")), 
             p1, 
             p2, 
             key, 
             textGrob("log(Cumulative Added Public Transit Trips)", gp=gpar(fontsize=20, font=2)),
             nrow=5, 
             heights=c(0.02, 0.42, 0.42, 0.08, 0.06)),
             rectGrob(gp=gpar(col="white")),
             arrangeGrob(us_cumul_plot,
                         arrangeGrob(arrangeGrob(tmax.hist, gcm_plot, ncol=2),
                                     arrangeGrob(rectGrob(gp=gpar(col="white")), 
                                                 justy.month, 
                                                 month20402050.plot, 
                                                 month209002100.plot, 
                                                 ncol=4, 
                                                 widths=c(0.035, 0.065, 0.45, 0.45)),
                                     nrow=2),
                         ncol=2),
             nrow=3, heights=c(0.49, 0.02, 0.49))
grid.text("a", x=unit(0.02, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30, font=2))
grid.text("b", x=unit(0.06, "npc"), y=unit(0.48, "npc"), gp=gpar(fontsize=30, font=2))
grid.text("c", x=unit(0.575, "npc"), y=unit(0.48, "npc"), gp=gpar(fontsize=30, font=2))
grid.text("d", x=unit(0.8175, "npc"), y=unit(0.48, "npc"), gp=gpar(fontsize=30, font=2))
grid.text("e", x=unit(0.57, "npc"), y=unit(0.06, "npc"), gp=gpar(fontsize=30, font=2))
grid.text("f", x=unit(0.795, "npc"), y=unit(0.06, "npc"), gp=gpar(fontsize=30, font=2))
grid.text("2050", x=unit(0.5, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30, font=2))
grid.text("2099", x=unit(0.5, "npc"), y=unit(0.775, "npc"), gp=gpar(fontsize=30, font=2))
grid.text("RCP2.6", x=unit(0.1675, "npc"), y=unit(0.9, "npc"), gp=gpar(fontsize=30, font=2, col="white"))
grid.text("RCP2.6", x=unit(0.1675, "npc"), y=unit(0.6925, "npc"), gp=gpar(fontsize=30, font=2, col="white"))
grid.text("RCP4.5", x=unit(0.4922, "npc"), y=unit(0.9, "npc"), gp=gpar(fontsize=30, font=2, col="white"))
grid.text("RCP4.5", x=unit(0.4922, "npc"), y=unit(0.6925, "npc"), gp=gpar(fontsize=30, font=2, col="white"))
grid.text("RCP8.5", x=unit(0.815, "npc"), y=unit(0.9, "npc"), gp=gpar(fontsize=30, font=2, col="white"))
grid.text("RCP8.5", x=unit(0.815, "npc"), y=unit(0.6925, "npc"), gp=gpar(fontsize=30, font=2, col="white"))
dev.off()
```

```{r ptt.projection.fig3, out.width="100%", fig.cap=paste0(sfig_num("ptt.fig3", display="cite"),": Climate change may amplify future public transit trips. Panel (a) plots the by-state log of projected cumulative added PTT by 2050 and 2099 across emissions scenarios. We map metropolitan areas into the states where the centroid of the metropolitan boundary area falls. Wyoming lacks baseline data during this period and so is omitted from the projection. Panel (b) displays the time-series of projected cumulative added PTT due to future warming over the entire US. Under RCP8.5, future US warming may produce nearly six billion cumulative added PTT by 2100. Panel (c) plots historical maximum monthly temperatures from our sample and the peak of the relationship between maximum temperatures and log(PTT). Approximately 90\\% of historical temperatures fall below this peak. Panel (d) plots monthly temperature anomaly projections relative to the 2006-2020 baseline for each state from 2020 to 2100. Panels (e) and (f) plot sumtotal added PTT during each month of the year for 2040-2050 and 2090-2100.")}
knitr::include_graphics("./ptt_projection.png")
```

```{r}
ptt.fig3 <- sfig_num(name = "ptt.fig3", caption = "")
```

```{r ptt.monthlyproj, eval=FALSE}
# calculate and plot 2040-2050 and 2090-2100 by-city map projections----
## note, collapse cities to states by the centroid of their metro area
uza_projection <- readRDS("./uza_projection_bootstrap.rds")

## need to merge in state ids based on uza centroids
names(states)[6] <- "state" # fix names
states$state <- as.character(states$state)
points <- uza_projection[, .(lat = unique(lat), lon = unique(lon)), by=.(uza)]
coordinates(points) <- ~lon + lat
projection(points) <- projection(states)
uza.state <- cbind(as.data.table(points), as.data.table(over(points, states)))
uza.state <- uza.state[, .(uza, state)]
setkey(uza.state, uza)
setkey(uza_projection, uza)
uza_projection <- merge(uza_projection, uza.state)
uza_projection[uza==13, state := "California"]
uza_projection[uza==17, state := "Florida"]
uza_projection[uza==34, state := "Virginia"]

## average uza_projection to state-level
uza_projection <- uza_projection[, lapply(.SD, mean), .SDcols=c("tmax_rcp45","tmax_rcp85","tmax_rcp26","tmax_rcp26_coef","tmax_rcp45_coef","tmax_rcp85_coef","tmax_rcp26_upr","tmax_rcp26_lwr","tmax_rcp45_upr","tmax_rcp45_lwr","tmax_rcp85_upr","tmax_rcp85_lwr"),
                                 by=.(state, monthyr, month)]

## need to use 2006-2020 tmax baseline and calculate differences by state then cumulate
mean2006_2020 <- uza_projection[year(monthyr) >= 2006 & year(monthyr) < 2020, 
                                  .(base_rcp26_coef = mean(tmax_rcp26_coef),
                                    base_rcp45_coef = mean(tmax_rcp45_coef),
                                    base_rcp85_coef = mean(tmax_rcp85_coef)), 
                                  by=.(month(monthyr), state)]
uza_projection[, month := month(monthyr)]
setkey(uza_projection, state, month)
setkey(mean2006_2020, state, month)
uza_projection <- merge(uza_projection, mean2006_2020)
rm(mean2006_2020)

## get baseline ptt by state for 2008-2017 (most recent full ten-year period)
setkey(uza.state, uza)
setkey(upt, uza)
upt <- merge(upt, uza.state, all.x=T)
baseline_ptt <- upt[, .(trips = sum(trips, na.rm=T)), 
                    by=.(state, monthyr)] # sum to state-level
baseline_ptt <- baseline_ptt[year(monthyr) >= 2008 & year(monthyr) <= 2017, 
                            .(trips = mean(trips, na.rm=T)), 
                            by=.(state, month = month(monthyr))]
setkey(baseline_ptt, state, month)
setkey(uza_projection, state, month)
uza_projection <- merge(uza_projection, baseline_ptt)
rm(baseline_ptt)

## calculate changes relative to baseline
## exponentiate coefficients to get percentage change, difference from fixed baseline
uza_projection[, ':=' (delta_rcp26 = 100*(exp(tmax_rcp26_coef) - exp(base_rcp26_coef)),
                         delta_rcp45 = 100*(exp(tmax_rcp45_coef) - exp(base_rcp45_coef)),
                         delta_rcp85 = 100*(exp(tmax_rcp85_coef) - exp(base_rcp85_coef)))] 

## keep only year 2020+
uza_projection <- uza_projection[year(monthyr) >= 2020]

## get by-month
month20402050 <- uza_projection[year(monthyr) >= 2040 & year(monthyr) <= 2050, 
                                 .(delta_rcp26 = mean(delta_rcp26, na.rm=T),
                                   delta_rcp45 = mean(delta_rcp45, na.rm=T),
                                   delta_rcp85 = mean(delta_rcp85, na.rm=T)),
                                 by=.(month, state=as.character(state))]
month20402050 <- melt(month20402050, id.vars = c('month', 'state'), variable.name = "rcp", value.name = 'trips')

month20902100 <- uza_projection[year(monthyr) >= 2090 & year(monthyr) <= 2100, 
                                 .(delta_rcp26 = mean(delta_rcp26, na.rm=T),
                                   delta_rcp45 = mean(delta_rcp45, na.rm=T),
                                   delta_rcp85 = mean(delta_rcp85, na.rm=T)),
                                 by=.(month, state=as.character(state))]
month20902100 <- melt(month20902100, id.vars = c('month','state'), variable.name = "rcp", value.name = 'trips')

## plot data
names(states)[6] <- "state" # fix names
states$state <- as.character(states$state)
month2050rcp26 <- month20402050[rcp=="delta_rcp26", .(trips), by=.(month, state)]
month2099rcp26 <- month20902100[rcp=="delta_rcp26", .(trips), by=.(month, state)]
month2050rcp26 <- dcast(month2050rcp26, state ~ paste0('month',month), value.var="trips")
month2099rcp26 <- dcast(month2099rcp26, state ~ paste0('month',month), value.var="trips")
month2050rcp26 <- sp::merge(states, month2050rcp26, by="state", duplicateGeoms=TRUE)
month2099rcp26 <- sp::merge(states, month2099rcp26, by="state", duplicateGeoms=TRUE)

month2050rcp45 <- month20402050[rcp=="delta_rcp45", .(trips), by=.(month, state)]
month2099rcp45 <- month20902100[rcp=="delta_rcp45", .(trips), by=.(month, state)]
month2050rcp45 <- dcast(month2050rcp45, state ~ paste0('month',month), value.var="trips")
month2099rcp45 <- dcast(month2099rcp45, state ~ paste0('month',month), value.var="trips")
month2050rcp45 <- sp::merge(states, month2050rcp45, by="state", duplicateGeoms=TRUE)
month2099rcp45 <- sp::merge(states, month2099rcp45, by="state", duplicateGeoms=TRUE)

month2050rcp85 <- month20402050[rcp=="delta_rcp85", .(trips), by=.(month, state)]
month2099rcp85 <- month20902100[rcp=="delta_rcp85", .(trips), by=.(month, state)]
month2050rcp85 <- dcast(month2050rcp85, state ~ paste0('month',month), value.var="trips")
month2099rcp85 <- dcast(month2099rcp85, state ~ paste0('month',month), value.var="trips")
month2050rcp85 <- sp::merge(states, month2050rcp85, by="state", duplicateGeoms=TRUE)
month2099rcp85 <- sp::merge(states, month2099rcp85, by="state", duplicateGeoms=TRUE)

# rcp26----
p1 <- spplot(
       obj=month2050rcp26,
       zcol=c('month1','month2','month3','month4','month5','month6','month7','month8','month9','month10','month11','month12'),         
       names.attr = paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December")), # set names for strip text
       strip=TRUE, # include strip text
       col.regions=colorRampPalette(c("#d5ed76","#8ed1ff","#56bbff","#19A1FC", "#0083db"))(100), # set colors
       as.table = TRUE,
       # col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2), height=0.5), # make legend on bottom of plot
       at=seq(min(as.data.table(cbind(month2050rcp26[11:ncol(month2050rcp26)], month2099rcp85[11:ncol(month2099rcp85)])), na.rm=T), max(as.data.table(cbind(month2050rcp26[11:ncol(month2050rcp26)], month2099rcp85[11:ncol(month2099rcp85)])), na.rm=T)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       layout=c(4,3) # horizontal layout
       )
p1 <- update(p1, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p1 <- p1 + layer_(sp.lines(obj=us, fill="white", col="black"))
p1 <- p1 + layer(sp.lines(obj=states,
                              fill="transparent",
                              col="black",
                              cex=0.1))

p2 <- spplot(
       obj=month2099rcp26,
       zcol=c('month1','month2','month3','month4','month5','month6','month7','month8','month9','month10','month11','month12'),         
       names.attr = paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December")), # set names for strip text
       strip=TRUE, # include strip text
       col.regions=colorRampPalette(c("#d5ed76","#8ed1ff","#56bbff","#19A1FC", "#0083db"))(100), # set colors
       as.table = TRUE,
       # col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2), height=0.5), # make legend on bottom of plot
       at=seq(min(as.data.table(cbind(month2050rcp26[11:ncol(month2050rcp26)], month2099rcp85[11:ncol(month2099rcp85)])), na.rm=T), max(as.data.table(cbind(month2050rcp26[11:ncol(month2050rcp26)], month2099rcp85[11:ncol(month2099rcp85)])), na.rm=T)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       layout=c(4,3) # horizontal layout
       )
p2 <- update(p2, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p2 <- p2 + layer_(sp.lines(obj=us, fill="white", col="black"))
p2 <- p2 + layer(sp.lines(obj=states,
                              fill="transparent",
                              col="black",
                              cex=0.1))

## get legend for plotting
key <- draw.colorkey(p1$legend[[1]]$args$key)
p1$legend <- NULL
p2$legend <- NULL
p1$sub <- NULL
p2$sub <- NULL

png(filename = './ptt_bymonth_rcp26.png', width=9.75, height=11, units = "in", type="cairo", res=300)
grid.arrange(rectGrob(gp=gpar(col="white")), 
             p1, 
             p2, 
             key, 
             textGrob("Percentage Change in PTT, RCP2.6", gp=gpar(fontsize=23, font=2)),
             nrow=5, 
             heights=c(0.02, 0.43, 0.43, 0.07, 0.05))
grid.text("2040-2050", x=unit(0.5, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=21, font=2))
grid.text("2090-2100", x=unit(0.5, "npc"), y=unit(0.545, "npc"), gp=gpar(fontsize=21, font=2))
dev.off()

# rcp45----
p1 <- spplot(
       obj=month2050rcp45,
       zcol=c('month1','month2','month3','month4','month5','month6','month7','month8','month9','month10','month11','month12'),         
       names.attr = paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December")), # set names for strip text
       strip=TRUE, # include strip text
       col.regions=colorRampPalette(c("#d5ed76","#8ed1ff","#56bbff","#19A1FC", "#0083db"))(100), # set colors
       as.table = TRUE,
       # col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2), height=0.5), # make legend on bottom of plot
       at=seq(min(as.data.table(cbind(month2050rcp26[11:ncol(month2050rcp26)], month2099rcp85[11:ncol(month2099rcp85)])), na.rm=T), max(as.data.table(cbind(month2050rcp26[11:ncol(month2050rcp26)], month2099rcp85[11:ncol(month2099rcp85)])), na.rm=T)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       layout=c(4,3) # horizontal layout
       )
p1 <- update(p1, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p1 <- p1 + layer_(sp.lines(obj=us, fill="white", col="black"))
p1 <- p1 + layer(sp.lines(obj=states,
                              fill="transparent",
                              col="black",
                              cex=0.1))

p2 <- spplot(
       obj=month2099rcp45,
       zcol=c('month1','month2','month3','month4','month5','month6','month7','month8','month9','month10','month11','month12'),         
       names.attr = paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December")), # set names for strip text
       strip=TRUE, # include strip text
       col.regions=colorRampPalette(c("#d5ed76","#8ed1ff","#56bbff","#19A1FC", "#0083db"))(100), # set colors
       as.table = TRUE,
       # col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2), height=0.5), # make legend on bottom of plot
       at=seq(min(as.data.table(cbind(month2050rcp26[11:ncol(month2050rcp26)], month2099rcp85[11:ncol(month2099rcp85)])), na.rm=T), max(as.data.table(cbind(month2050rcp26[11:ncol(month2050rcp26)], month2099rcp85[11:ncol(month2099rcp85)])), na.rm=T)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       layout=c(4,3) # horizontal layout
       )
p2 <- update(p2, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p2 <- p2 + layer_(sp.lines(obj=us, fill="white", col="black"))
p2 <- p2 + layer(sp.lines(obj=states,
                              fill="transparent",
                              col="black",
                              cex=0.1))

## get legend for plotting
key <- draw.colorkey(p1$legend[[1]]$args$key)
p1$legend <- NULL
p2$legend <- NULL
p1$sub <- NULL
p2$sub <- NULL

png(filename = './ptt_bymonth_rcp45.png', width=9.75, height=11, units = "in", type="cairo", res=300)
grid.arrange(rectGrob(gp=gpar(col="white")), 
             p1, 
             p2, 
             key, 
             textGrob("Percentage Change in PTT, RCP4.5", gp=gpar(fontsize=23, font=2)),
             nrow=5, 
             heights=c(0.02, 0.43, 0.43, 0.07, 0.05))
grid.text("2040-2050", x=unit(0.5, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=21, font=2))
grid.text("2090-2100", x=unit(0.5, "npc"), y=unit(0.545, "npc"), gp=gpar(fontsize=21, font=2))
dev.off()

# rcp85----
p1 <- spplot(
       obj=month2050rcp85,
       zcol=c('month1','month2','month3','month4','month5','month6','month7','month8','month9','month10','month11','month12'),         
       names.attr = paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December")), # set names for strip text
       strip=TRUE, # include strip text
       col.regions=colorRampPalette(c("#d5ed76","#8ed1ff","#56bbff","#19A1FC", "#0083db"))(100), # set colors
       as.table = TRUE,
       # col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2), height=0.5), # make legend on bottom of plot
       at=seq(min(as.data.table(cbind(month2050rcp26[11:ncol(month2050rcp26)], month2099rcp85[11:ncol(month2099rcp85)])), na.rm=T), max(as.data.table(cbind(month2050rcp26[11:ncol(month2050rcp26)], month2099rcp85[11:ncol(month2099rcp85)])), na.rm=T)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       layout=c(4,3) # horizontal layout
       )
p1 <- update(p1, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p1 <- p1 + layer_(sp.lines(obj=us, fill="white", col="black"))
p1 <- p1 + layer(sp.lines(obj=states,
                              fill="transparent",
                              col="black",
                              cex=0.1))

p2 <- spplot(
       obj=month2099rcp85,
       zcol=c('month1','month2','month3','month4','month5','month6','month7','month8','month9','month10','month11','month12'),         
       names.attr = paste(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October","November","December")), # set names for strip text
       strip=TRUE, # include strip text
       col.regions=colorRampPalette(c("#d5ed76","#8ed1ff","#56bbff","#19A1FC", "#0083db"))(100), # set colors
       as.table = TRUE,
       # col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2), height=0.5), # make legend on bottom of plot
       at=seq(min(as.data.table(cbind(month2050rcp26[11:ncol(month2050rcp26)], month2099rcp85[11:ncol(month2099rcp85)])), na.rm=T), max(as.data.table(cbind(month2050rcp26[11:ncol(month2050rcp26)], month2099rcp85[11:ncol(month2099rcp85)])), na.rm=T)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       layout=c(4,3) # horizontal layout
       )
p2 <- update(p2, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p2 <- p2 + layer_(sp.lines(obj=us, fill="white", col="black"))
p2 <- p2 + layer(sp.lines(obj=states,
                              fill="transparent",
                              col="black",
                              cex=0.1))

## get legend for plotting
key <- draw.colorkey(p1$legend[[1]]$args$key)
p1$legend <- NULL
p2$legend <- NULL
p1$sub <- NULL
p2$sub <- NULL

png(filename = './ptt_bymonth_rcp85.png', width=9.75, height=11, units = "in", type="cairo", res=300)
grid.arrange(rectGrob(gp=gpar(col="white")), 
             p1, 
             p2, 
             key, 
             textGrob("Percentage Change in PTT, RCP8.5", gp=gpar(fontsize=23, font=2)),
             nrow=5, 
             heights=c(0.02, 0.43, 0.43, 0.07, 0.05))
grid.text("2040-2050", x=unit(0.5, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=21, font=2))
grid.text("2090-2100", x=unit(0.5, "npc"), y=unit(0.545, "npc"), gp=gpar(fontsize=21, font=2))
dev.off()
```

```{r ptt.monthlyproj26, out.width="100%", fig.cap = paste0(sfig_num("ptt.monthlyproj26", display="cite"),": Plot of monthly projected mean percentage changes in public transit trips 2040-2050 and 2090-2100 across the RCP2.6 emissions scenario. This figure replicates the methods of Figure 4 from the main text across the RCP2.6 emissions scenario for PTT. We map metropolitan areas into the states where the centroid of the metropolitan boundary area falls. Wyoming lacks baseline data during this period and so is omitted from the projection. The RCP2.6 emissions scenario results in lower projected changes in future PTT.")}
knitr::include_graphics("./ptt_bymonth_rcp26.png")
```

```{r}
ptt.monthlyproj26 <- sfig_num(name = "ptt.monthlyproj26", caption = "")
```

```{r ptt.monthlyproj45, out.width="100%", fig.cap = paste0(sfig_num("ptt.monthlyproj45", display="cite"),": Plot of monthly projected mean percentage changes in public transit trips 2040-2050 and 2090-2100 across the RCP4.5 emissions scenario. This figure replicates the methods of Figure 4 from the main text across the RCP4.5 emissions scenario for PTT. We map metropolitan areas into the states where the centroid of the metropolitan boundary area falls. Wyoming lacks baseline data during this period and so is omitted from the projection. The RCP4.5 emissions scenario results in middle-range projected changes in future PTT.")}
knitr::include_graphics("./ptt_bymonth_rcp45.png")
```

```{r}
ptt.monthlyproj45 <- sfig_num(name = "ptt.monthlyproj45", caption = "")
```

```{r ptt.monthlyproj85, out.width="100%", fig.cap = paste0(sfig_num("ptt.monthlyproj85", display="cite"),": Plot of monthly projected mean percentage changes in public transit trips 2040-2050 and 2090-2100 across the RCP8.5 emissions scenario. This figure replicates the methods of Figure 4 from the main text across the RCP8.5 emissions scenario for PTT. We map metropolitan areas into the states where the centroid of the metropolitan boundary area falls. Wyoming lacks baseline data during this period and so is omitted from the projection. The RCP8.5 emissions scenario results in higher-range projected changes in future PTT.")}
knitr::include_graphics("./ptt_bymonth_rcp85.png")
```

```{r}
ptt.monthlyproj85 <- sfig_num(name = "ptt.monthlyproj85", caption = "")
```

<!-- tables -->

```{r spline.reg}
# VMT----
data <- vmt
dv <- vmt$log_vmt
iv <- vmt$tmax
lhs <- "log_vmt"

## keep spline knots with the fitted model results for plotting
spl_knots <- wtd.quantile(iv[!is.na(dv)], weights=dv[!is.na(dv)], probs = c(0.25, 0.5, 0.75), na.rm=T) # get quantiles of iv weighted by dv (cannot have missing values)
spl_boundary_knots <- range(iv, na.rm=T) # get range of independent variable
xs <- as.data.table(ns(iv, knots = spl_knots, Boundary.knots = spl_boundary_knots)) # get cubic splines of iv
setnames(xs, paste0("ns_", 1:ncol(xs))) # set names of splines
data <- cbind(data, xs)
## set up felm model
rhs <- paste(c(names(xs), "trange", "prcp.days", "cloud", "humid", "wind"), collapse = " + ")
fe <- paste(c("state:month", "state:year"), collapse = " + ")
fmla <- as.formula(sprintf("%s ~ %s | %s | 0 | state", lhs, rhs, fe))
r1 <- felm(fmla, data=data, exactDOF = TRUE)

# PTT----
data <- upt
dv <- upt$log_trips
iv <- upt$tmax
lhs <- "log_trips"

## keep spline knots with the fitted model results for plotting
spl_knots <- wtd.quantile(iv[!is.na(dv)], weights=dv[!is.na(dv)], probs = c(0.25, 0.5, 0.75), na.rm=T) # get quantiles of iv weighted by dv (cannot have missing values)
spl_boundary_knots <- range(iv, na.rm=T) # get range of independent variable
xs <- as.data.table(ns(iv, knots = spl_knots, Boundary.knots = spl_boundary_knots)) # get cubic splines of iv
setnames(xs, paste0("ns_", 1:ncol(xs))) # set names of splines
data <- cbind(data, xs)
## set up felm model
rhs <- paste(c(names(xs), "trange", "prcp.days", "cloud", "humid", "wind"), collapse = " + ")
fe <- paste(c("uza:month", "uza:year"), collapse = " + ")
fmla <- as.formula(sprintf("%s ~ %s | %s | 0 | uza", lhs, rhs, fe))
r2 <- felm(fmla, data=data, exactDOF=TRUE)
```

```{r spline.tab, results='asis'}
invisible(stargazer(r1, r2,
          type='latex',
          title = paste0(stab_num("spline.reg", display="cite"), ':', " Splined Regressions from Equation 1 for VMT and PTT"),
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c("VMT", "PTT"),
          covariate.labels = c("TMAX SPLINE 1", "TMAX SPLINE 2", "TMAX SPLINE 3", "TMAX SPLINE 4", "TRANGE", "PRCP.DAYS","HUMID","CLOUD","WIND"),
          add.lines = list(c("State:Month FE", "Yes", "No"),
                           c("State:Year FE", "Yes", "No"), 
                           c("City:Month FE", "No", "Yes"),
                           c("City:Year FE", "No", "Yes")),
          notes = c('Standard errors in parentheses are clustered on state (VMT) and city (PTT).'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          table.placement = "h"
          ))
```

```{r}
spline.reg <- stab_num(name = "spline.reg", caption = "")
```
